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, (char *)0, 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(VecSet(vecLocal1, 0.0));
35: PetscCall(DMGlobalToLocalBegin(dm, vec, INSERT_VALUES, vecLocal1));
36: PetscCall(DMGlobalToLocalEnd(dm, vec, INSERT_VALUES, vecLocal1));
38: PetscCall(DMStagGetCorners(dm, &startx, &starty, &startz, &nx, &ny, &nz, NULL, NULL, NULL));
39: PetscCall(DMStagVecGetArrayRead(dm, vecLocal1, &a1));
40: PetscCall(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) PetscCall(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: PetscCall(DMStagVecRestoreArrayRead(dm, vecLocal1, &a1));
56: PetscCall(DMStagVecRestoreArray(dm, vecLocal2, &a2));
58: PetscCall(DMLocalToGlobalBegin(dm, vecLocal2, INSERT_VALUES, vec));
59: PetscCall(DMLocalToGlobalEnd(dm, vecLocal2, INSERT_VALUES, vec));
61: /* For the all-periodic case, some additional checks */
62: PetscCall(DMStagGetBoundaryTypes(dm, &boundaryTypex, &boundaryTypey, &boundaryTypez));
63: if (boundaryTypex == DM_BOUNDARY_PERIODIC && boundaryTypey == DM_BOUNDARY_PERIODIC && boundaryTypez == DM_BOUNDARY_PERIODIC) {
64: PetscCall(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: PetscCall(VecSum(vecLocal1, &sum));
67: 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)));
69: PetscCall(VecGetArray(vec, &a));
70: expected = 1 + 6 * stencilWidth;
71: for (i = 0; i < nz * ny * nx * dofTotal; ++i) {
72: if (a[i] != expected) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Unexpected value %g (expecting %g)\n", rank, (double)PetscRealPart(a[i]), (double)PetscRealPart(expected)));
73: }
74: PetscCall(VecRestoreArray(vec, &a));
75: }
77: PetscCall(VecDestroy(&vec));
78: PetscCall(VecDestroy(&vecLocal1));
79: PetscCall(VecDestroy(&vecLocal2));
80: PetscCall(DMDestroy(&dm));
81: PetscCall(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*/