xref: /petsc/src/vec/vec/tutorials/ex9.c (revision 15229ffc342989b2bf0590a733d92c152a3348fc)
1 static char help[] = "Demonstrates use of VecCreateGhost().\n\n";
2 
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 */
10 
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>
18 
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;
28 
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");
34 
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
39 
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().
44 
45      In this example we consider the vector as representing
46      degrees of freedom in a one dimensional grid with periodic
47      boundary conditions.
48 
49         ----Processor  1---------  ----Processor 2 --------
50          0    1   2   3   4    5    6    7   8   9   10   11
51                                |----|
52          |-------------------------------------------------|
53 
54   */
55 
56   if (rank == 0) {
57     ifrom[0] = 11;
58     ifrom[1] = 6;
59   } else {
60     ifrom[0] = 0;
61     ifrom[1] = 5;
62   }
63 
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   }
83 
84   /*
85       Test VecDuplicate()
86   */
87   PetscCall(VecDuplicate(gxs, &gx));
88   PetscCall(VecDestroy(&gxs));
89 
90   /*
91      Access the local representation
92   */
93   PetscCall(VecGhostGetLocalForm(gx, &lx));
94 
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));
105 
106   PetscCall(VecGhostUpdateBegin(gx, INSERT_VALUES, SCATTER_FORWARD));
107   PetscCall(VecGhostUpdateEnd(gx, INSERT_VALUES, SCATTER_FORWARD));
108 
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));
117 
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));
126 
127     PetscCall(VecGhostUpdateBegin(gx, MIN_VALUES, SCATTER_REVERSE));
128     PetscCall(VecGhostUpdateEnd(gx, MIN_VALUES, SCATTER_REVERSE));
129 
130     PetscCall(VecGhostGetLocalForm(gx, &lx));
131     PetscCall(VecGetArray(lx, &array));
132 
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   }
138 
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   }
149 
150   PetscCall(VecDestroy(&gx));
151 
152   if (flg) PetscCall(PetscFree(tarray));
153   PetscCall(PetscFinalize());
154   return 0;
155 }
156 
157 /*TEST
158 
159      test:
160        nsize: 2
161 
162      test:
163        suffix: 2
164        nsize: 2
165        args: -allocate
166        output_file: output/ex9_1.out
167 
168      test:
169        suffix: 3
170        nsize: 2
171        args: -vecmpisetghost
172        output_file: output/ex9_1.out
173 
174      test:
175        suffix: 4
176        nsize: 2
177        args: -minvalues
178        output_file: output/ex9_2.out
179        requires: !complex
180 
181      test:
182        suffix: 5
183        nsize: 2
184        args: -vecghostgetghostis
185 
186      test:
187        suffix: 6
188        nsize: 2
189        args: -getgtlmapping
190 
191 TEST*/
192