Actual source code: ex9.c

  1: static char help[] = "Test DMStag 3d star stencil\n\n";
  2: #include <petscdm.h>
  3: #include <petscdmstag.h>

  5: int main(int argc, char **argv)
  6: {
  7:   DM             dm;
  8:   Vec            vec, vecLocal1, vecLocal2;
  9:   PetscScalar   *a, ****a1, ****a2, expected, sum;
 10:   PetscInt       startx, starty, startz, nx, ny, nz, i, j, k, d, is, js, ks, dof0, dof1, dof2, dof3, dofTotal, stencilWidth, ngx, ngy, ngz;
 11:   DMBoundaryType boundaryTypex, boundaryTypey, boundaryTypez;
 12:   PetscMPIInt    rank;

 15:   PetscInitialize(&argc, &argv, (char *)0, help);
 16:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 17:   dof0         = 1;
 18:   dof1         = 1;
 19:   dof2         = 1;
 20:   dof3         = 1;
 21:   stencilWidth = 2;
 22:   DMStagCreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, 4, 4, 4, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, dof0, dof1, dof2, dof3, DMSTAG_STENCIL_STAR, stencilWidth, NULL, NULL, NULL, &dm);
 23:   DMSetFromOptions(dm);
 24:   DMSetUp(dm);
 25:   DMStagGetDOF(dm, &dof0, &dof1, &dof2, &dof3);
 26:   dofTotal = dof0 + 3 * dof1 + 3 * dof2 + dof3;
 27:   DMStagGetStencilWidth(dm, &stencilWidth);

 29:   DMCreateLocalVector(dm, &vecLocal1);
 30:   VecDuplicate(vecLocal1, &vecLocal2);

 32:   DMCreateGlobalVector(dm, &vec);
 33:   VecSet(vec, 1.0);
 34:   VecSet(vecLocal1, 0.0);
 35:   DMGlobalToLocalBegin(dm, vec, INSERT_VALUES, vecLocal1);
 36:   DMGlobalToLocalEnd(dm, vec, INSERT_VALUES, vecLocal1);

 38:   DMStagGetCorners(dm, &startx, &starty, &startz, &nx, &ny, &nz, NULL, NULL, NULL);
 39:   DMStagVecGetArrayRead(dm, vecLocal1, &a1);
 40:   DMStagVecGetArray(dm, vecLocal2, &a2);
 41:   for (k = startz; k < startz + nz; ++k) {
 42:     for (j = starty; j < starty + ny; ++j) {
 43:       for (i = startx; i < startx + nx; ++i) {
 44:         for (d = 0; d < dofTotal; ++d) {
 45:           if (a1[k][j][i][d] != 1.0) PetscPrintf(PETSC_COMM_SELF, "[%d] Unexpected value %g (expecting %g)\n", rank, (double)PetscRealPart(a1[k][j][i][d]), 1.0);
 46:           a2[k][j][i][d] = 0.0;
 47:           for (ks = -stencilWidth; ks <= stencilWidth; ++ks) a2[k][j][i][d] += a1[k + ks][j][i][d];
 48:           for (js = -stencilWidth; js <= stencilWidth; ++js) a2[k][j][i][d] += a1[k][j + js][i][d];
 49:           for (is = -stencilWidth; is <= stencilWidth; ++is) a2[k][j][i][d] += a1[k][j][i + is][d];
 50:           a2[k][j][i][d] -= 2.0 * a1[k][j][i][d];
 51:         }
 52:       }
 53:     }
 54:   }
 55:   DMStagVecRestoreArrayRead(dm, vecLocal1, &a1);
 56:   DMStagVecRestoreArray(dm, vecLocal2, &a2);

 58:   DMLocalToGlobalBegin(dm, vecLocal2, INSERT_VALUES, vec);
 59:   DMLocalToGlobalEnd(dm, vecLocal2, INSERT_VALUES, vec);

 61:   /* For the all-periodic case, some additional checks */
 62:   DMStagGetBoundaryTypes(dm, &boundaryTypex, &boundaryTypey, &boundaryTypez);
 63:   if (boundaryTypex == DM_BOUNDARY_PERIODIC && boundaryTypey == DM_BOUNDARY_PERIODIC && boundaryTypez == DM_BOUNDARY_PERIODIC) {
 64:     DMStagGetGhostCorners(dm, NULL, NULL, NULL, &ngx, &ngy, &ngz);
 65:     expected = (ngx * ngy * ngz - 8 * stencilWidth * stencilWidth * stencilWidth - 4 * stencilWidth * stencilWidth * (nx + ny + nz)) * dofTotal;
 66:     VecSum(vecLocal1, &sum);
 67:     if (sum != expected) PetscPrintf(PETSC_COMM_SELF, "[%d] Unexpected sum of local entries %g (expected %g)\n", rank, (double)PetscRealPart(sum), (double)PetscRealPart(expected));

 69:     VecGetArray(vec, &a);
 70:     expected = 1 + 6 * stencilWidth;
 71:     for (i = 0; i < nz * ny * nx * dofTotal; ++i) {
 72:       if (a[i] != expected) PetscPrintf(PETSC_COMM_SELF, "[%d] Unexpected value %g (expecting %g)\n", rank, (double)PetscRealPart(a[i]), (double)PetscRealPart(expected));
 73:     }
 74:     VecRestoreArray(vec, &a);
 75:   }

 77:   VecDestroy(&vec);
 78:   VecDestroy(&vecLocal1);
 79:   VecDestroy(&vecLocal2);
 80:   DMDestroy(&dm);
 81:   PetscFinalize();
 82:   return 0;
 83: }

 85: /*TEST

 87:    test:
 88:       suffix: 1
 89:       nsize: 8
 90:       args: -stag_ranks_x 2 -stag_ranks_y 2 -stag_ranks_z 2 -stag_stencil_width 1

 92:    test:
 93:       suffix: 2
 94:       nsize: 12
 95:       args: -stag_ranks_x 3 -stag_ranks_y 2 -stag_ranks_z 2 -stag_dof_0 2 -stag_grid_x 6

 97:    test:
 98:       suffix: 3
 99:       nsize: 8
100:       args: -stag_dof_0 3 -stag_dof_1 2 -stag_dof_2 4 -stag_dof_3 2 -stag_stencil_width 3 -stag_grid_x 6 -stag_grid_y 6 -stag_grid_z 7

102:    test:
103:       suffix: 4
104:       nsize: 8
105:       args: -stag_stencil_width 1 -stag_grid_x 2 -stag_grid_y 2 -stag_grid_z 2 -stag_boundary_type_x ghosted

107:    test:
108:       suffix: 5
109:       nsize: 8
110:       args: -stag_stencil_width 1 -stag_grid_x 2 -stag_grid_y 2 -stag_grid_z 2 -stag_boundary_type_y ghosted

112:    test:
113:       suffix: 6
114:       nsize: 8
115:       args: -stag_stencil_width 1 -stag_grid_x 2 -stag_grid_y 2 -stag_grid_z 2 -stag_boundary_type_z ghosted -stag_dof_0 2 -stag_dof_1 3 -stag_dof_2 2 -stag_dof_3 2

117:    test:
118:       suffix: 7
119:       nsize: 8
120:       args: -stag_stencil_width 1 -stag_grid_x 3 -stag_grid_y 2 -stag_grid_z 2 -stag_boundary_type_x ghosted -stag_boundary_type_y ghosted

122:    test:
123:       suffix: 8
124:       nsize: 8
125:       args: -stag_stencil_width 1 -stag_grid_x 2 -stag_grid_y 5 -stag_grid_z 2 -stag_boundary_type_x ghosted -stag_boundary_type_z ghosted

127:    test:
128:       suffix: 9
129:       nsize: 8
130:       args: -stag_stencil_width 1 -stag_grid_x 2 -stag_grid_y 2 -stag_grid_z 3 -stag_boundary_type_y ghosted -stag_boundary_type_z ghosted

132:    test:
133:       suffix: 10
134:       nsize: 5
135:       args: -stag_stencil_width 1 -stag_grid_y 2 -stag_grid_z 3 -stag_grid_x 17 -stag_boundary_type_y ghosted -stag_boundary_type_z ghosted -stag_ranks_x 5
136: TEST*/