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*/