1 static char help[] = "Test DMStag 3d star stencil\n\n";
2 #include <petscdm.h>
3 #include <petscdmstag.h>
4
main(int argc,char ** argv)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;
13
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));
28
29 PetscCall(DMCreateLocalVector(dm, &vecLocal1));
30 PetscCall(VecDuplicate(vecLocal1, &vecLocal2));
31
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));
37
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));
57
58 PetscCall(DMLocalToGlobalBegin(dm, vecLocal2, INSERT_VALUES, vec));
59 PetscCall(DMLocalToGlobalEnd(dm, vecLocal2, INSERT_VALUES, vec));
60
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)));
68
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 }
76
77 PetscCall(VecDestroy(&vec));
78 PetscCall(VecDestroy(&vecLocal1));
79 PetscCall(VecDestroy(&vecLocal2));
80 PetscCall(DMDestroy(&dm));
81 PetscCall(PetscFinalize());
82 return 0;
83 }
84
85 /*TEST
86
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
91 output_file: output/empty.out
92
93 test:
94 suffix: 2
95 nsize: 12
96 args: -stag_ranks_x 3 -stag_ranks_y 2 -stag_ranks_z 2 -stag_dof_0 2 -stag_grid_x 6
97 output_file: output/empty.out
98
99 test:
100 suffix: 3
101 nsize: 8
102 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
103 output_file: output/empty.out
104
105 test:
106 suffix: 4
107 nsize: 8
108 args: -stag_stencil_width 1 -stag_grid_x 2 -stag_grid_y 2 -stag_grid_z 2 -stag_boundary_type_x ghosted
109 output_file: output/empty.out
110
111 test:
112 suffix: 5
113 nsize: 8
114 args: -stag_stencil_width 1 -stag_grid_x 2 -stag_grid_y 2 -stag_grid_z 2 -stag_boundary_type_y ghosted
115 output_file: output/empty.out
116
117 test:
118 suffix: 6
119 nsize: 8
120 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
121 output_file: output/empty.out
122
123 test:
124 suffix: 7
125 nsize: 8
126 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
127 output_file: output/empty.out
128
129 test:
130 suffix: 8
131 nsize: 8
132 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
133 output_file: output/empty.out
134
135 test:
136 suffix: 9
137 nsize: 8
138 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
139 output_file: output/empty.out
140
141 test:
142 suffix: 10
143 nsize: 5
144 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
145 output_file: output/empty.out
146 TEST*/
147