Actual source code: ex49.c

petsc-3.11.3 2019-06-26
Report Typos and Errors
  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*/