Actual source code: ex49.c

  1: static const char help[] = "Test VEC_SUBSET_OFF_PROC_ENTRIES\n\n";

  3: #include <petsc.h>
  4: #include <petscvec.h>

  6: /* Unlike most finite element applications, IBAMR does assembly on many cells
  7:    that are not locally owned; in some cases the processor may own zero finite
  8:    element cells but still do assembly on a small number of cells anyway. To
  9:    simulate this, this code assembles a PETSc vector by adding contributions
 10:    to every entry in the vector on every processor. This causes a deadlock
 11:    when we save the communication pattern via

 13:      VecSetOption(vec, VEC_SUBSET_OFF_PROC_ENTRIES, PETSC_TRUE).

 15:    Contributed-by: David Wells <drwells@email.unc.edu>

 17:   Petsc developers' notes: this test tests how Petsc knows it can reuse existing communication
 18:   pattern. All processes must come to the same conclusion, otherwise deadlock may happen due
 19:   to mismatched MPI_Send/Recv. It also tests changing VEC_SUBSET_OFF_PROC_ENTRIES back and forth.
 20: */
 21: int main(int argc, char **argv)
 22: {
 23:   Vec         v;
 24:   PetscInt    i, j, k, *ln, n, rstart;
 25:   PetscBool   saveCommunicationPattern = PETSC_FALSE;
 26:   PetscMPIInt size, rank, p;

 29:   PetscInitialize(&argc, &argv, NULL, help);
 30:   PetscOptionsGetBool(NULL, NULL, "-save_comm", &saveCommunicationPattern, NULL);
 31:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 32:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

 34:   PetscMalloc1(size, &ln);
 35:   /* This bug is triggered when one of the local lengths is small. Sometimes in IBAMR this value is actually zero. */
 36:   for (p = 0; p < size; ++p) ln[p] = 10;
 37:   ln[0] = 2;
 38:   PetscPrintf(PETSC_COMM_WORLD, "local lengths are:\n");
 39:   PetscIntView(1, &ln[rank], PETSC_VIEWER_STDOUT_WORLD);
 40:   n = ln[rank];
 41:   VecCreateMPI(MPI_COMM_WORLD, n, PETSC_DECIDE, &v);
 42:   VecGetOwnershipRange(v, &rstart, NULL);

 44:   for (k = 0; k < 5; ++k) { /* 5 iterations of VecAssembly */
 45:     PetscReal norm  = 0.0;
 46:     PetscBool flag  = (k == 2) ? PETSC_FALSE : PETSC_TRUE;
 47:     PetscInt  shift = (k < 2) ? 0 : (k == 2) ? 1 : 0; /* Used to change patterns */

 49:     /* If saveCommunicationPattern, let's see what should happen in the 5 iterations:
 50:       iter 0: flag is true, and this is the first assembly, so petsc should keep the
 51:               communication pattern built during this assembly.
 52:       iter 1: flag is true, reuse the pattern.
 53:       iter 2: flag is false, discard/free the pattern built in iter 0; rebuild a new
 54:               pattern, but do not keep it after VecAssemblyEnd since the flag is false.
 55:       iter 3: flag is true again, this is the new first assembly with a true flag. So
 56:               petsc should keep the communication pattern built during this assembly.
 57:       iter 4: flag is true, reuse the pattern built in iter 3.

 59:       When the vector is destroyed, memory used by the pattern is freed. One can also do it early with a call
 60:           VecSetOption(v, VEC_SUBSET_OFF_PROC_ENTRIES, PETSC_FALSE);
 61:      */
 62:     if (saveCommunicationPattern) VecSetOption(v, VEC_SUBSET_OFF_PROC_ENTRIES, flag);
 63:     VecSet(v, 0.0);

 65:     for (i = 0; i < n; ++i) {
 66:       PetscScalar val = 1.0;
 67:       PetscInt    r   = rstart + i;

 69:       VecSetValue(v, r, val, ADD_VALUES);
 70:       /* do assembly on all other processors too (the 'neighbors') */
 71:       {
 72:         const PetscMPIInt neighbor = (i + shift) % size; /* Adjust communication patterns between iterations */
 73:         const PetscInt    nn       = ln[neighbor];
 74:         PetscInt          nrstart  = 0;

 76:         for (p = 0; p < neighbor; ++p) nrstart += ln[p];
 77:         for (j = 0; j < nn / 4; j += 3) {
 78:           PetscScalar val = 0.01;
 79:           PetscInt    nr  = nrstart + j;

 81:           VecSetValue(v, nr, val, ADD_VALUES);
 82:         }
 83:       }
 84:     }
 85:     VecAssemblyBegin(v);
 86:     VecAssemblyEnd(v);
 87:     VecNorm(v, NORM_1, &norm);
 88:     PetscPrintf(PETSC_COMM_WORLD, "norm is %g\n", (double)norm);
 89:   }
 90:   PetscFree(ln);
 91:   VecDestroy(&v);
 92:   PetscFinalize();
 93:   return 0;
 94: }

 96: /*TEST

 98:    test:
 99:       suffix: 1
100:       nsize: 4
101:    test:
102:       suffix: 1_save
103:       args: -save_comm
104:       nsize: 4
105:       output_file: output/ex49_1.out
106: TEST*/