Actual source code: ex4.c


  2: static char help[] = "Tests various 2-dimensional DMDA routines.\n\n";

  4: #include <petscdm.h>
  5: #include <petscdmda.h>

  7: int main(int argc, char **argv)
  8: {
  9:   PetscMPIInt     rank;
 10:   PetscInt        M = 10, N = 8, m = PETSC_DECIDE;
 11:   PetscInt        s = 2, w = 2, n = PETSC_DECIDE, nloc, l, i, j, kk;
 12:   PetscInt        Xs, Xm, Ys, Ym, iloc, *iglobal;
 13:   const PetscInt *ltog;
 14:   PetscInt       *lx = NULL, *ly = NULL;
 15:   PetscBool       testorder = PETSC_FALSE, flg;
 16:   DMBoundaryType  bx = DM_BOUNDARY_NONE, by = DM_BOUNDARY_NONE;
 17:   DM              da;
 18:   PetscViewer     viewer;
 19:   Vec             local, global;
 20:   PetscScalar     value;
 21:   DMDAStencilType st = DMDA_STENCIL_BOX;
 22:   AO              ao;

 25:   PetscInitialize(&argc, &argv, (char *)0, help);
 26:   PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, "", 300, 0, 400, 400, &viewer);

 28:   /* Readoptions */
 29:   PetscOptionsGetInt(NULL, NULL, "-NX", &M, NULL);
 30:   PetscOptionsGetInt(NULL, NULL, "-NY", &N, NULL);
 31:   PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
 32:   PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
 33:   PetscOptionsGetInt(NULL, NULL, "-s", &s, NULL);
 34:   PetscOptionsGetInt(NULL, NULL, "-w", &w, NULL);

 36:   flg = PETSC_FALSE;
 37:   PetscOptionsGetBool(NULL, NULL, "-xperiodic", &flg, NULL);
 38:   if (flg) bx = DM_BOUNDARY_PERIODIC;
 39:   flg = PETSC_FALSE;
 40:   PetscOptionsGetBool(NULL, NULL, "-yperiodic", &flg, NULL);
 41:   if (flg) by = DM_BOUNDARY_PERIODIC;
 42:   flg = PETSC_FALSE;
 43:   PetscOptionsGetBool(NULL, NULL, "-xghosted", &flg, NULL);
 44:   if (flg) bx = DM_BOUNDARY_GHOSTED;
 45:   flg = PETSC_FALSE;
 46:   PetscOptionsGetBool(NULL, NULL, "-yghosted", &flg, NULL);
 47:   if (flg) by = DM_BOUNDARY_GHOSTED;
 48:   flg = PETSC_FALSE;
 49:   PetscOptionsGetBool(NULL, NULL, "-star", &flg, NULL);
 50:   if (flg) st = DMDA_STENCIL_STAR;
 51:   flg = PETSC_FALSE;
 52:   PetscOptionsGetBool(NULL, NULL, "-box", &flg, NULL);
 53:   if (flg) st = DMDA_STENCIL_BOX;
 54:   flg = PETSC_FALSE;
 55:   PetscOptionsGetBool(NULL, NULL, "-testorder", &testorder, NULL);
 56:   /*
 57:       Test putting two nodes in x and y on each processor, exact last processor
 58:       in x and y gets the rest.
 59:   */
 60:   flg = PETSC_FALSE;
 61:   PetscOptionsGetBool(NULL, NULL, "-distribute", &flg, NULL);
 62:   if (flg) {
 64:     PetscMalloc1(m, &lx);
 65:     for (i = 0; i < m - 1; i++) lx[i] = 4;
 66:     lx[m - 1] = M - 4 * (m - 1);
 68:     PetscMalloc1(n, &ly);
 69:     for (i = 0; i < n - 1; i++) ly[i] = 2;
 70:     ly[n - 1] = N - 2 * (n - 1);
 71:   }

 73:   /* Create distributed array and get vectors */
 74:   DMDACreate2d(PETSC_COMM_WORLD, bx, by, st, M, N, m, n, w, s, lx, ly, &da);
 75:   DMSetFromOptions(da);
 76:   DMSetUp(da);
 77:   PetscFree(lx);
 78:   PetscFree(ly);

 80:   DMView(da, viewer);
 81:   DMCreateGlobalVector(da, &global);
 82:   DMCreateLocalVector(da, &local);

 84:   /* Set global vector; send ghost points to local vectors */
 85:   value = 1;
 86:   VecSet(global, value);
 87:   DMGlobalToLocalBegin(da, global, INSERT_VALUES, local);
 88:   DMGlobalToLocalEnd(da, global, INSERT_VALUES, local);

 90:   /* Scale local vectors according to processor rank; pass to global vector */
 91:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 92:   value = rank;
 93:   VecScale(local, value);
 94:   DMLocalToGlobalBegin(da, local, INSERT_VALUES, global);
 95:   DMLocalToGlobalEnd(da, local, INSERT_VALUES, global);

 97:   if (!testorder) { /* turn off printing when testing ordering mappings */
 98:     PetscPrintf(PETSC_COMM_WORLD, "\nGlobal Vectors:\n");
 99:     VecView(global, PETSC_VIEWER_STDOUT_WORLD);
100:     PetscPrintf(PETSC_COMM_WORLD, "\n\n");
101:   }

103:   /* Send ghost points to local vectors */
104:   DMGlobalToLocalBegin(da, global, INSERT_VALUES, local);
105:   DMGlobalToLocalEnd(da, global, INSERT_VALUES, local);

107:   flg = PETSC_FALSE;
108:   PetscOptionsGetBool(NULL, NULL, "-local_print", &flg, NULL);
109:   if (flg) {
110:     PetscViewer sviewer;

112:     PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);
113:     PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\nLocal Vector: processor %d\n", rank);
114:     PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer);
115:     VecView(local, sviewer);
116:     PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer);
117:     PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);
118:     PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD);
119:   }

121:   /* Tests mappings between application/PETSc orderings */
122:   if (testorder) {
123:     ISLocalToGlobalMapping ltogm;

125:     DMGetLocalToGlobalMapping(da, &ltogm);
126:     ISLocalToGlobalMappingGetSize(ltogm, &nloc);
127:     ISLocalToGlobalMappingGetIndices(ltogm, &ltog);
128:     DMDAGetGhostCorners(da, &Xs, &Ys, NULL, &Xm, &Ym, NULL);
129:     DMDAGetAO(da, &ao);
130:     PetscMalloc1(nloc, &iglobal);

132:     /* Set iglobal to be global indices for each processor's local and ghost nodes,
133:        using the DMDA ordering of grid points */
134:     kk = 0;
135:     for (j = Ys; j < Ys + Ym; j++) {
136:       for (i = Xs; i < Xs + Xm; i++) {
137:         iloc = w * ((j - Ys) * Xm + i - Xs);
138:         for (l = 0; l < w; l++) iglobal[kk++] = ltog[iloc + l];
139:       }
140:     }

142:     /* Map this to the application ordering (which for DMDAs is just the natural ordering
143:        that would be used for 1 processor, numbering most rapidly by x, then y) */
144:     AOPetscToApplication(ao, nloc, iglobal);

146:     /* Then map the application ordering back to the PETSc DMDA ordering */
147:     AOApplicationToPetsc(ao, nloc, iglobal);

149:     /* Verify the mappings */
150:     kk = 0;
151:     for (j = Ys; j < Ys + Ym; j++) {
152:       for (i = Xs; i < Xs + Xm; i++) {
153:         iloc = w * ((j - Ys) * Xm + i - Xs);
154:         for (l = 0; l < w; l++) {
155:           if (iglobal[kk] != ltog[iloc + l]) {
156:             PetscFPrintf(PETSC_COMM_SELF, stdout, "[%d] Problem with mapping: j=%" PetscInt_FMT ", i=%" PetscInt_FMT ", l=%" PetscInt_FMT ", petsc1=%" PetscInt_FMT ", petsc2=%" PetscInt_FMT "\n", rank, j, i, l, ltog[iloc + l], iglobal[kk]);
157:           }
158:           kk++;
159:         }
160:       }
161:     }
162:     PetscFree(iglobal);
163:     ISLocalToGlobalMappingRestoreIndices(ltogm, &ltog);
164:   }

166:   /* Free memory */
167:   PetscViewerDestroy(&viewer);
168:   VecDestroy(&local);
169:   VecDestroy(&global);
170:   DMDestroy(&da);

172:   PetscFinalize();
173:   return 0;
174: }

176: /*TEST

178:    test:
179:       nsize: 4
180:       args: -nox
181:       filter: grep -v -i Object
182:       requires: x

184:    test:
185:       suffix: 2
186:       args: -testorder -nox
187:       requires: x

189: TEST*/