1 static char help[] = "Test DMStag 2d 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, nx, ny, i, j, d, is, js, dof0, dof1, dof2, dofTotal, stencilWidth, ngx, ngy;
11 DMBoundaryType boundaryTypex, boundaryTypey;
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 stencilWidth = 2;
21 PetscCall(DMStagCreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, 4, 4, PETSC_DECIDE, PETSC_DECIDE, dof0, dof1, dof2, DMSTAG_STENCIL_STAR, stencilWidth, NULL, NULL, &dm));
22 PetscCall(DMSetFromOptions(dm));
23 PetscCall(DMSetUp(dm));
24 PetscCall(DMStagGetDOF(dm, &dof0, &dof1, &dof2, NULL));
25 dofTotal = dof0 + 2 * dof1 + dof2;
26 PetscCall(DMStagGetStencilWidth(dm, &stencilWidth));
27
28 PetscCall(DMCreateLocalVector(dm, &vecLocal1));
29 PetscCall(VecDuplicate(vecLocal1, &vecLocal2));
30
31 PetscCall(DMCreateGlobalVector(dm, &vec));
32 PetscCall(VecSet(vec, 1.0));
33 PetscCall(VecSet(vecLocal1, 0.0));
34 PetscCall(DMGlobalToLocalBegin(dm, vec, INSERT_VALUES, vecLocal1));
35 PetscCall(DMGlobalToLocalEnd(dm, vec, INSERT_VALUES, vecLocal1));
36
37 PetscCall(DMStagGetCorners(dm, &startx, &starty, NULL, &nx, &ny, NULL, NULL, NULL, NULL));
38 PetscCall(DMStagVecGetArrayRead(dm, vecLocal1, &a1));
39 PetscCall(DMStagVecGetArray(dm, vecLocal2, &a2));
40 for (j = starty; j < starty + ny; ++j) {
41 for (i = startx; i < startx + nx; ++i) {
42 for (d = 0; d < dofTotal; ++d) {
43 if (a1[j][i][d] != 1.0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Unexpected value %g (expecting %g)\n", rank, (double)PetscRealPart(a1[j][i][d]), 1.0));
44 a2[j][i][d] = 0.0;
45 for (js = -stencilWidth; js <= stencilWidth; ++js) a2[j][i][d] += a1[j + js][i][d];
46 for (is = -stencilWidth; is <= stencilWidth; ++is) a2[j][i][d] += a1[j][i + is][d];
47 a2[j][i][d] -= a1[j][i][d];
48 }
49 }
50 }
51 PetscCall(DMStagVecRestoreArrayRead(dm, vecLocal1, &a1));
52 PetscCall(DMStagVecRestoreArray(dm, vecLocal2, &a2));
53
54 PetscCall(DMLocalToGlobalBegin(dm, vecLocal2, INSERT_VALUES, vec));
55 PetscCall(DMLocalToGlobalEnd(dm, vecLocal2, INSERT_VALUES, vec));
56
57 /* For the all-periodic case, some additional checks */
58 PetscCall(DMStagGetBoundaryTypes(dm, &boundaryTypex, &boundaryTypey, NULL));
59 if (boundaryTypex == DM_BOUNDARY_PERIODIC && boundaryTypey == DM_BOUNDARY_PERIODIC) {
60 PetscCall(DMStagGetGhostCorners(dm, NULL, NULL, NULL, &ngx, &ngy, NULL));
61 expected = (ngx * ngy - 4 * stencilWidth * stencilWidth) * dofTotal;
62 PetscCall(VecSum(vecLocal1, &sum));
63 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)));
64
65 PetscCall(VecGetArray(vec, &a));
66 expected = 1 + 4 * stencilWidth;
67 for (i = 0; i < ny * nx * dofTotal; ++i) {
68 if (a[i] != expected) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Unexpected value %g (expecting %g)\n", rank, (double)PetscRealPart(a[i]), (double)PetscRealPart(expected)));
69 }
70 PetscCall(VecRestoreArray(vec, &a));
71 }
72
73 PetscCall(VecDestroy(&vec));
74 PetscCall(VecDestroy(&vecLocal1));
75 PetscCall(VecDestroy(&vecLocal2));
76 PetscCall(DMDestroy(&dm));
77 PetscCall(PetscFinalize());
78 return 0;
79 }
80
81 /*TEST
82
83 test:
84 suffix: 1
85 nsize: 4
86 args: -stag_ranks_x 2 -stag_ranks_y 2 -stag_stencil_width 1
87 output_file: output/empty.out
88
89 test:
90 suffix: 2
91 nsize: 6
92 args: -stag_ranks_x 3 -stag_ranks_y 2 -stag_dof_0 2 -stag_grid_x 6
93 output_file: output/empty.out
94
95 test:
96 suffix: 3
97 nsize: 4
98 args: -stag_dof_0 3 -stag_dof_1 2 -stag_dof_2 4 2 -stag_stencil_width 3 -stag_grid_x 6 -stag_grid_y 6
99 output_file: output/empty.out
100
101 test:
102 suffix: 4
103 nsize: 4
104 args: -stag_stencil_width 1 -stag_grid_x 2 -stag_grid_y 2 -stag_boundary_type_x ghosted
105 output_file: output/empty.out
106
107 test:
108 suffix: 5
109 nsize: 4
110 args: -stag_stencil_width 1 -stag_grid_x 2 -stag_grid_y 2 -stag_boundary_type_y ghosted
111 output_file: output/empty.out
112
113 test:
114 suffix: 6
115 nsize: 4
116 args: -stag_stencil_width 1 -stag_grid_x 3 -stag_grid_y 2 -stag_boundary_type_x ghosted -stag_boundary_type_y ghosted
117 output_file: output/empty.out
118
119 test:
120 suffix: 7
121 nsize: 4
122 args: -stag_stencil_width 1 -stag_grid_x 2 -stag_grid_y 2 -stag_boundary_type_y ghosted
123 output_file: output/empty.out
124
125 test:
126 suffix: 8
127 nsize: 6
128 args: -stag_stencil_width 1 -stag_grid_y 2 -stag_grid_x 19 -stag_boundary_type_y ghosted -stag_ranks_x 6
129 output_file: output/empty.out
130 TEST*/
131