Actual source code: ex49.c
petsc-3.11.4 2019-09-28
1: static const char help[] = "Test VEC_SUBSET_OFF_PROC_ENTRIES\n\n";
3: #include <petsc.h>
4: #include <petscvec.h>
6: #include <stdlib.h>
8: /* Unlike most finite element applications, IBAMR does assembly on many cells
9: that are not locally owned; in some cases the processor may own zero finite
10: element cells but still do assembly on a small number of cells anyway. To
11: simulate this, this code assembles a PETSc vector by adding contributions
12: to every entry in the vector on every processor. This causes a deadlock
13: when we save the communication pattern via
15: VecSetOption(vec, VEC_SUBSET_OFF_PROC_ENTRIES, PETSC_TRUE).
17: Contributed-by: David Wells <drwells@email.unc.edu>
19: Petsc developers' notes: this test tests how Petsc knows it can reuse existing communication
20: pattern. All processes must come to the same conclusion, otherwise deadlock may happen due
21: to mismatched MPI_Send/Recv. It also tests changing VEC_SUBSET_OFF_PROC_ENTRIES back and forth.
22: */
23: int main(int argc, char **argv)
24: {
25: Vec v;
26: PetscInt i, j, k, *ln, n, rstart;
27: PetscBool saveCommunicationPattern = PETSC_FALSE;
28: PetscMPIInt size, rank, p;
31: PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
32: PetscOptionsGetBool(NULL, NULL, "-save_comm", &saveCommunicationPattern, NULL);
33: MPI_Comm_size(PETSC_COMM_WORLD, &size);
34: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
36: PetscMalloc1(size, &ln);
37: /* This bug is triggered when one of the local lengths is small. Sometimes in IBAMR this value is actually zero. */
38: for (p=0; p<size; ++p) ln[p] = 10;
39: ln[0] = 2;
40: PetscPrintf(PETSC_COMM_WORLD, "local lengths are:\n");
41: PetscIntView(1, &ln[rank], PETSC_VIEWER_STDOUT_WORLD);
42: n = ln[rank];
43: VecCreateMPI(MPI_COMM_WORLD, n, PETSC_DECIDE, &v);
44: VecGetOwnershipRange(v, &rstart, NULL);
46: for (k=0; k<5; ++k) { /* 5 iterations of VecAssembly */
47: PetscReal norm = 0.0;
48: PetscBool flag = (k == 2) ? PETSC_FALSE : PETSC_TRUE;
49: PetscInt shift = (k < 2) ? 0 : (k == 2) ? 1 : 0; /* Used to change patterns */
51: /* If saveCommunicationPattern, let's see what should happen in the 5 iterations:
52: iter 0: flag is true, and this is the first assebmly, so petsc should keep the
53: communication pattern built during this assembly.
54: iter 1: flag is true, reuse the pattern.
55: iter 2: flag is false, discard/free the pattern built in iter 0; rebuild a new
56: pattern, but do not keep it after VecAssemblyEnd since the flag is false.
57: iter 3: flag is true again, this is the new first assembly with a true flag. So
58: petsc should keep the communication pattern built during this assembly.
59: iter 4: flag is true, reuse the pattern built in iter 3.
61: When the vector is destroyed, memory used by the pattern is freed. One can also do it early with a call
62: VecSetOption(v, VEC_SUBSET_OFF_PROC_ENTRIES, PETSC_FALSE);
63: */
64: if (saveCommunicationPattern) {VecSetOption(v, VEC_SUBSET_OFF_PROC_ENTRIES, flag);}
65: VecSet(v, 0.0);
67: for (i=0; i<n; ++i) {
68: PetscScalar val = 1.0;
69: PetscInt r = rstart + i;
71: VecSetValue(v, r, val, ADD_VALUES);
72: /* do assembly on all other processors too (the 'neighbors') */
73: {
74: const PetscMPIInt neighbor = (i+shift) % size; /* Adjust communication patterns between iterations */
75: const PetscInt nn = ln[neighbor];
76: PetscInt nrstart = 0;
78: for (p=0; p<neighbor; ++p) nrstart += ln[p];
79: for (j=0; j<nn/4; j+= 3) {
80: PetscScalar val = 0.01;
81: PetscInt nr = nrstart + j;
83: VecSetValue(v, nr, val, ADD_VALUES);
84: }
85: }
86: }
87: VecAssemblyBegin(v);
88: VecAssemblyEnd(v);
89: VecNorm(v, NORM_1, &norm);
90: PetscPrintf(PETSC_COMM_WORLD, "norm is %g\n", (double)norm);
91: }
92: PetscFree(ln);
93: VecDestroy(&v);
94: PetscFinalize();
95: return ierr;
96: }
98: /*TEST
100: test:
101: suffix: 1
102: nsize: 4
103: test:
104: suffix: 1_save
105: args: -save_comm
106: nsize: 4
107: output_file: output/ex49_1.out
108: TEST*/