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;
14: PetscFunctionBeginUser;
15: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
16: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
17: dof0 = 1;
18: dof1 = 1;
19: dof2 = 1;
20: dof3 = 1;
21: stencilWidth = 2;
22: PetscCall(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: PetscCall(DMSetFromOptions(dm));
24: PetscCall(DMSetUp(dm));
25: PetscCall(DMStagGetDOF(dm, &dof0, &dof1, &dof2, &dof3));
26: dofTotal = dof0 + 3 * dof1 + 3 * dof2 + dof3;
27: PetscCall(DMStagGetStencilWidth(dm, &stencilWidth));
29: PetscCall(DMCreateLocalVector(dm, &vecLocal1));
30: PetscCall(VecDuplicate(vecLocal1, &vecLocal2));
32: PetscCall(DMCreateGlobalVector(dm, &vec));
33: PetscCall(VecSet(vec, 1.0));
34: PetscCall(DMGlobalToLocalBegin(dm, vec, INSERT_VALUES, vecLocal1));
35: PetscCall(DMGlobalToLocalEnd(dm, vec, INSERT_VALUES, vecLocal1));
37: PetscCall(DMStagGetCorners(dm, &startx, &starty, &startz, &nx, &ny, &nz, NULL, NULL, NULL));
38: PetscCall(DMStagVecGetArrayRead(dm, vecLocal1, &a1));
39: PetscCall(DMStagVecGetArray(dm, vecLocal2, &a2));
40: for (k = startz; k < startz + nz; ++k) {
41: for (j = starty; j < starty + ny; ++j) {
42: for (i = startx; i < startx + nx; ++i) {
43: for (d = 0; d < dofTotal; ++d) {
44: if (a1[k][j][i][d] != 1.0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Unexpected value %g (expecting %g)\n", rank, (double)PetscRealPart(a1[k][j][i][d]), 1.0));
45: a2[k][j][i][d] = 0.0;
46: for (ks = -stencilWidth; ks <= stencilWidth; ++ks) a2[k][j][i][d] += a1[k + ks][j][i][d];
47: for (js = -stencilWidth; js <= stencilWidth; ++js) a2[k][j][i][d] += a1[k][j + js][i][d];
48: for (is = -stencilWidth; is <= stencilWidth; ++is) a2[k][j][i][d] += a1[k][j][i + is][d];
49: a2[k][j][i][d] -= 2.0 * a1[k][j][i][d];
50: }
51: }
52: }
53: }
54: PetscCall(DMStagVecRestoreArrayRead(dm, vecLocal1, &a1));
55: PetscCall(DMStagVecRestoreArray(dm, vecLocal2, &a2));
57: PetscCall(DMLocalToGlobalBegin(dm, vecLocal2, INSERT_VALUES, vec));
58: PetscCall(DMLocalToGlobalEnd(dm, vecLocal2, INSERT_VALUES, vec));
60: /* For the all-periodic case, some additional checks */
61: PetscCall(DMStagGetBoundaryTypes(dm, &boundaryTypex, &boundaryTypey, &boundaryTypez));
62: if (boundaryTypex == DM_BOUNDARY_PERIODIC && boundaryTypey == DM_BOUNDARY_PERIODIC && boundaryTypez == DM_BOUNDARY_PERIODIC) {
63: PetscCall(DMStagGetGhostCorners(dm, NULL, NULL, NULL, &ngx, &ngy, &ngz));
64: expected = (ngx * ngy * ngz - 8 * stencilWidth * stencilWidth * stencilWidth - 4 * stencilWidth * stencilWidth * (nx + ny + nz)) * dofTotal;
65: PetscCall(VecSum(vecLocal1, &sum));
66: if (sum != expected) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Unexpected sum of local entries %g (expected %g)\n", rank, (double)PetscRealPart(sum), (double)PetscRealPart(expected)));
68: PetscCall(VecGetArray(vec, &a));
69: expected = 1 + 6 * stencilWidth;
70: for (i = 0; i < nz * ny * nx * dofTotal; ++i) {
71: if (a[i] != expected) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Unexpected value %g (expecting %g)\n", rank, (double)PetscRealPart(a[i]), (double)PetscRealPart(expected)));
72: }
73: PetscCall(VecRestoreArray(vec, &a));
74: }
76: PetscCall(VecDestroy(&vec));
77: PetscCall(VecDestroy(&vecLocal1));
78: PetscCall(VecDestroy(&vecLocal2));
79: PetscCall(DMDestroy(&dm));
80: PetscCall(PetscFinalize());
81: return 0;
82: }
84: /*TEST
86: test:
87: suffix: 1
88: nsize: 8
89: args: -stag_ranks_x 2 -stag_ranks_y 2 -stag_ranks_z 2 -stag_stencil_width 1
90: output_file: output/empty.out
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
96: output_file: output/empty.out
98: test:
99: suffix: 3
100: nsize: 8
101: 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: output_file: output/empty.out
104: test:
105: suffix: 4
106: nsize: 8
107: args: -stag_stencil_width 1 -stag_grid_x 2 -stag_grid_y 2 -stag_grid_z 2 -stag_boundary_type_x ghosted
108: output_file: output/empty.out
110: test:
111: suffix: 5
112: nsize: 8
113: args: -stag_stencil_width 1 -stag_grid_x 2 -stag_grid_y 2 -stag_grid_z 2 -stag_boundary_type_y ghosted
114: output_file: output/empty.out
116: test:
117: suffix: 6
118: nsize: 8
119: 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
120: output_file: output/empty.out
122: test:
123: suffix: 7
124: nsize: 8
125: 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
126: output_file: output/empty.out
128: test:
129: suffix: 8
130: nsize: 8
131: 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
132: output_file: output/empty.out
134: test:
135: suffix: 9
136: nsize: 8
137: 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
138: output_file: output/empty.out
140: test:
141: suffix: 10
142: nsize: 5
143: 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
144: output_file: output/empty.out
145: TEST*/