1 static char help[] = "Test DMStag 2d periodic and ghosted boundary conditions\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;
10 PetscInt startx, starty, nx, ny, i, j, d, is, js, dof0, dof1, dof2, dofTotal, stencilWidth, Nx, Ny;
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_BOX, 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) {
46 for (is = -stencilWidth; is <= stencilWidth; ++is) a2[j][i][d] += a1[j + js][i + is][d];
47 }
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, all values are the same . Otherwise, just check the local version */
58 PetscCall(DMStagGetBoundaryTypes(dm, &boundaryTypex, &boundaryTypey, NULL));
59 if (boundaryTypex == DM_BOUNDARY_PERIODIC && boundaryTypey == DM_BOUNDARY_PERIODIC) {
60 PetscCall(VecGetArray(vec, &a));
61 expected = 1.0;
62 for (d = 0; d < 2; ++d) expected *= (2 * stencilWidth + 1);
63 for (i = 0; i < ny * nx * dofTotal; ++i) {
64 if (a[i] != expected) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Unexpected value %g (expecting %g)\n", rank, (double)PetscRealPart(a[i]), (double)PetscRealPart(expected)));
65 }
66 PetscCall(VecRestoreArray(vec, &a));
67 } else {
68 PetscCall(DMStagVecGetArrayRead(dm, vecLocal2, &a2));
69 PetscCall(DMStagGetGlobalSizes(dm, &Nx, &Ny, NULL));
70 PetscCheck(stencilWidth <= 1, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Non-periodic check implemented assuming stencilWidth = 1");
71 for (j = starty; j < starty + ny; ++j) {
72 for (i = startx; i < startx + nx; ++i) {
73 PetscInt dd, extra[2];
74 PetscBool bnd[2];
75 bnd[0] = (PetscBool)((i == 0 || i == Nx - 1) && boundaryTypex != DM_BOUNDARY_PERIODIC);
76 bnd[1] = (PetscBool)((j == 0 || j == Ny - 1) && boundaryTypey != DM_BOUNDARY_PERIODIC);
77 extra[0] = i == Nx - 1 && boundaryTypex != DM_BOUNDARY_PERIODIC ? 1 : 0;
78 extra[1] = j == Ny - 1 && boundaryTypey != DM_BOUNDARY_PERIODIC ? 1 : 0;
79 { /* vertices */
80 PetscScalar expected = 1.0;
81 for (dd = 0; dd < 2; ++dd) expected *= (bnd[dd] ? stencilWidth + 1 + extra[dd] : 2 * stencilWidth + 1);
82 for (d = 0; d < dof0; ++d) {
83 if (a2[j][i][d] != expected) {
84 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Element (%" PetscInt_FMT ",%" PetscInt_FMT ")[%" PetscInt_FMT "] Unexpected value %g (expecting %g)\n", rank, i, j, d, (double)PetscRealPart(a2[j][i][d]), (double)PetscRealPart(expected)));
85 }
86 }
87 }
88 { /* down edges */
89 PetscScalar expected = (bnd[1] ? stencilWidth + 1 + extra[1] : 2 * stencilWidth + 1);
90 expected *= ((bnd[0] ? 1 : 2) * stencilWidth + 1);
91 for (d = dof0; d < dof0 + dof1; ++d) {
92 if (a2[j][i][d] != expected) {
93 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Element (%" PetscInt_FMT ",%" PetscInt_FMT ")[%" PetscInt_FMT "] Unexpected value %g (expecting %g)\n", rank, i, j, d, (double)PetscRealPart(a2[j][i][d]), (double)PetscRealPart(expected)));
94 }
95 }
96 }
97 { /* left edges */
98 PetscScalar expected = (bnd[0] ? stencilWidth + 1 + extra[0] : 2 * stencilWidth + 1);
99 expected *= ((bnd[1] ? 1 : 2) * stencilWidth + 1);
100 for (d = dof0 + dof1; d < dof0 + 2 * dof1; ++d) {
101 if (a2[j][i][d] != expected) {
102 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Element (%" PetscInt_FMT ",%" PetscInt_FMT ")[%" PetscInt_FMT "] Unexpected value %g (expecting %g)\n", rank, i, j, d, (double)PetscRealPart(a2[j][i][d]), (double)PetscRealPart(expected)));
103 }
104 }
105 }
106 { /* elements */
107 PetscScalar expected = 1.0;
108 for (dd = 0; dd < 2; ++dd) expected *= ((bnd[dd] ? 1 : 2) * stencilWidth + 1);
109 for (d = dofTotal - dof2; d < dofTotal; ++d) {
110 if (a2[j][i][d] != expected) {
111 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Element (%" PetscInt_FMT ",%" PetscInt_FMT ")[%" PetscInt_FMT "] Unexpected value %g (expecting %g)\n", rank, i, j, d, (double)PetscRealPart(a2[j][i][d]), (double)PetscRealPart(expected)));
112 }
113 }
114 }
115 }
116 }
117 PetscCall(DMStagVecRestoreArrayRead(dm, vecLocal2, &a2));
118 }
119
120 PetscCall(VecDestroy(&vec));
121 PetscCall(VecDestroy(&vecLocal1));
122 PetscCall(VecDestroy(&vecLocal2));
123 PetscCall(DMDestroy(&dm));
124 PetscCall(PetscFinalize());
125 return 0;
126 }
127
128 /*TEST
129
130 test:
131 suffix: 1
132 nsize: 4
133 args: -stag_ranks_x 2 -stag_ranks_y 2 -stag_stencil_width 1 -stag_dof_2 2
134 output_file: output/empty.out
135
136 test:
137 suffix: 2
138 nsize: 4
139 args: -stag_ranks_x 2 -stag_ranks_y 2 -stag_dof_1 2 -stag_grid_y 5
140 output_file: output/empty.out
141
142 test:
143 suffix: 3
144 nsize: 6
145 args: -stag_ranks_x 3 -stag_ranks_y 2 -stag_dof_0 2 -stag_grid_x 6
146 output_file: output/empty.out
147
148 test:
149 suffix: 4
150 nsize: 6
151 args: -stag_ranks_x 3 -stag_ranks_y 2 -stag_dof_0 0 -stag_dof_1 0 -stag_dof_2 0 -stag_grid_x 4 -stag_boundary_type_x ghosted -stag_boundary_type_y ghosted -stag_stencil_width 1
152 output_file: output/empty.out
153
154 test:
155 suffix: 5
156 nsize: 6
157 args: -stag_ranks_x 3 -stag_ranks_y 2 -stag_dof_0 0 -stag_dof_1 0 -stag_dof_2 0 -stag_grid_x 4 -stag_boundary_type_x ghosted -stag_stencil_width 1
158 output_file: output/empty.out
159
160 test:
161 suffix: 6
162 nsize: 9
163 args: -stag_dof_0 2 -stag_dof_1 2 -stag_dof_2 1 -stag_dof_2 1 -stag_boundary_type_y ghosted -stag_grid_x 9 -stag_grid_y 13 -stag_ranks_x 3 -stag_ranks_y 3 -stag_stencil_width 1
164 output_file: output/empty.out
165
166 test:
167 suffix: 7
168 nsize: 1
169 args: -stag_dof_0 2 -stag_dof_1 2 -stag_dof_2 1 -stag_dof_2 1 stag_grid_x 9 -stag_grid_y 13 -stag_stencil_width 1
170 output_file: output/empty.out
171
172 TEST*/
173