1 static char help[] = "Test DMStag 3d 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, startz, nx, ny, nz, i, j, k, d, is, js, ks, dof0, dof1, dof2, dof3, dofTotal, stencilWidth, Nx, Ny, Nz;
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_BOX, 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) {
48 for (js = -stencilWidth; js <= stencilWidth; ++js) {
49 for (is = -stencilWidth; is <= stencilWidth; ++is) a2[k][j][i][d] += a1[k + ks][j + js][i + is][d];
50 }
51 }
52 }
53 }
54 }
55 }
56 PetscCall(DMStagVecRestoreArrayRead(dm, vecLocal1, &a1));
57 PetscCall(DMStagVecRestoreArray(dm, vecLocal2, &a2));
58
59 PetscCall(DMLocalToGlobalBegin(dm, vecLocal2, INSERT_VALUES, vec));
60 PetscCall(DMLocalToGlobalEnd(dm, vecLocal2, INSERT_VALUES, vec));
61
62 /* For the all-periodic case, all values are the same . Otherwise, just check the local version */
63 PetscCall(DMStagGetBoundaryTypes(dm, &boundaryTypex, &boundaryTypey, &boundaryTypez));
64 if (boundaryTypex == DM_BOUNDARY_PERIODIC && boundaryTypey == DM_BOUNDARY_PERIODIC && boundaryTypez == DM_BOUNDARY_PERIODIC) {
65 PetscCall(VecGetArray(vec, &a));
66 expected = 1.0;
67 for (d = 0; d < 3; ++d) expected *= (2 * stencilWidth + 1);
68 for (i = 0; i < nz * ny * nx * dofTotal; ++i) {
69 if (a[i] != expected) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Unexpected value %g (expecting %g)\n", rank, (double)PetscRealPart(a[i]), (double)PetscRealPart(expected)));
70 }
71 PetscCall(VecRestoreArray(vec, &a));
72 } else {
73 PetscCall(DMStagVecGetArrayRead(dm, vecLocal2, &a2));
74 PetscCall(DMStagGetGlobalSizes(dm, &Nx, &Ny, &Nz));
75 PetscCheck(stencilWidth <= 1, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Check implemented assuming stencilWidth = 1");
76 for (k = startz; k < startz + nz; ++k) {
77 for (j = starty; j < starty + ny; ++j) {
78 for (i = startx; i < startx + nx; ++i) {
79 PetscInt dd, extra[3];
80 PetscBool bnd[3];
81 bnd[0] = (PetscBool)((i == 0 || i == Nx - 1) && boundaryTypex != DM_BOUNDARY_PERIODIC);
82 bnd[1] = (PetscBool)((j == 0 || j == Ny - 1) && boundaryTypey != DM_BOUNDARY_PERIODIC);
83 bnd[2] = (PetscBool)((k == 0 || k == Nz - 1) && boundaryTypez != DM_BOUNDARY_PERIODIC);
84 extra[0] = i == Nx - 1 && boundaryTypex != DM_BOUNDARY_PERIODIC ? 1 : 0;
85 extra[1] = j == Ny - 1 && boundaryTypey != DM_BOUNDARY_PERIODIC ? 1 : 0;
86 extra[2] = k == Nz - 1 && boundaryTypez != DM_BOUNDARY_PERIODIC ? 1 : 0;
87 { /* vertices */
88 PetscScalar expected = 1.0;
89 for (dd = 0; dd < 3; ++dd) expected *= (bnd[dd] ? stencilWidth + 1 + extra[dd] : 2 * stencilWidth + 1);
90 for (d = 0; d < dof0; ++d) {
91 if (a2[k][j][i][d] != expected) {
92 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Element (%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ")[%" PetscInt_FMT "] Unexpected value %g (expecting %g)\n", rank, i, j, k, d, (double)PetscRealPart(a2[k][j][i][d]), (double)PetscRealPart(expected)));
93 }
94 }
95 }
96 { /* back down edges */
97 PetscScalar expected = ((bnd[0] ? 1 : 2) * stencilWidth + 1);
98 for (dd = 1; dd < 3; ++dd) expected *= (bnd[dd] ? stencilWidth + 1 + extra[dd] : 2 * stencilWidth + 1);
99 for (d = dof0; d < dof0 + dof1; ++d) {
100 if (a2[k][j][i][d] != expected) {
101 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Element (%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ")[%" PetscInt_FMT "] Unexpected value %g (expecting %g)\n", rank, i, j, k, d, (double)PetscRealPart(a2[k][j][i][d]), (double)PetscRealPart(expected)));
102 }
103 }
104 }
105 { /* back left edges */
106 PetscScalar expected = ((bnd[1] ? 1 : 2) * stencilWidth + 1);
107 for (dd = 0; dd < 3; dd += 2) expected *= (bnd[dd] ? stencilWidth + 1 + extra[dd] : 2 * stencilWidth + 1);
108 for (d = dof0 + dof1; d < dof0 + 2 * dof1; ++d) {
109 if (a2[k][j][i][d] != expected) {
110 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Element (%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ")[%" PetscInt_FMT "] Unexpected value %g (expecting %g)\n", rank, i, j, k, d, (double)PetscRealPart(a2[k][j][i][d]), (double)PetscRealPart(expected)));
111 }
112 }
113 }
114 { /* back faces */
115 PetscScalar expected = (bnd[2] ? stencilWidth + 1 + extra[2] : 2 * stencilWidth + 1);
116 for (dd = 0; dd < 2; ++dd) expected *= ((bnd[dd] ? 1 : 2) * stencilWidth + 1);
117 for (d = dof0 + 2 * dof1; d < dof0 + 2 * dof1 + dof2; ++d) {
118 if (a2[k][j][i][d] != expected) {
119 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Element (%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ")[%" PetscInt_FMT "] Unexpected value %g (expecting %g)\n", rank, i, j, k, d, (double)PetscRealPart(a2[k][j][i][d]), (double)PetscRealPart(expected)));
120 }
121 }
122 }
123 { /* down left edges */
124 PetscScalar expected = ((bnd[2] ? 1 : 2) * stencilWidth + 1);
125 for (dd = 0; dd < 2; ++dd) expected *= (bnd[dd] ? stencilWidth + 1 + extra[dd] : 2 * stencilWidth + 1);
126 for (d = dof0 + 2 * dof1 + dof2; d < dof0 + 3 * dof1 + dof2; ++d) {
127 if (a2[k][j][i][d] != expected) {
128 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Element (%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ")[%" PetscInt_FMT "] Unexpected value %g (expecting %g)\n", rank, i, j, k, d, (double)PetscRealPart(a2[k][j][i][d]), (double)PetscRealPart(expected)));
129 }
130 }
131 }
132 { /* down faces */
133 PetscScalar expected = (bnd[1] ? stencilWidth + 1 + extra[1] : 2 * stencilWidth + 1);
134 for (dd = 0; dd < 3; dd += 2) expected *= ((bnd[dd] ? 1 : 2) * stencilWidth + 1);
135 for (d = dof0 + 3 * dof1 + dof2; d < dof0 + 3 * dof1 + 2 * dof2; ++d) {
136 if (a2[k][j][i][d] != expected) {
137 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Element (%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ")[%" PetscInt_FMT "] Unexpected value %g (expecting %g)\n", rank, i, j, k, d, (double)PetscRealPart(a2[k][j][i][d]), (double)PetscRealPart(expected)));
138 }
139 }
140 }
141 { /* left faces */
142 PetscScalar expected = (bnd[0] ? stencilWidth + 1 + extra[0] : 2 * stencilWidth + 1);
143 for (dd = 1; dd < 3; ++dd) expected *= ((bnd[dd] ? 1 : 2) * stencilWidth + 1);
144 for (d = dof0 + 3 * dof1 + 2 * dof2; d < dof0 + 3 * dof1 + 3 * dof2; ++d) {
145 if (a2[k][j][i][d] != expected) {
146 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Element (%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ")[%" PetscInt_FMT "] Unexpected value %g (expecting %g)\n", rank, i, j, k, d, (double)PetscRealPart(a2[k][j][i][d]), (double)PetscRealPart(expected)));
147 }
148 }
149 }
150 { /* elements */
151 PetscScalar expected = 1.0;
152 for (dd = 0; dd < 3; ++dd) expected *= ((bnd[dd] ? 1 : 2) * stencilWidth + 1);
153 for (d = dofTotal - dof3; d < dofTotal; ++d) {
154 if (a2[k][j][i][d] != expected) {
155 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Element (%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ")[%" PetscInt_FMT "] Unexpected value %g (expecting %g)\n", rank, i, j, k, d, (double)PetscRealPart(a2[k][j][i][d]), (double)PetscRealPart(expected)));
156 }
157 }
158 }
159 }
160 }
161 }
162 PetscCall(DMStagVecRestoreArrayRead(dm, vecLocal2, &a2));
163 }
164
165 PetscCall(VecDestroy(&vec));
166 PetscCall(VecDestroy(&vecLocal1));
167 PetscCall(VecDestroy(&vecLocal2));
168 PetscCall(DMDestroy(&dm));
169 PetscCall(PetscFinalize());
170 return 0;
171 }
172
173 /*TEST
174
175 test:
176 suffix: 1
177 nsize: 8
178 args: -stag_ranks_x 2 -stag_ranks_y 2 -stag_ranks_z 2 -stag_stencil_width 1 -stag_dof_3 2 -stag_grid_z 3
179 output_file: output/empty.out
180
181 test:
182 suffix: 2
183 nsize: 8
184 args: -stag_ranks_x 2 -stag_ranks_y 2 -stag_ranks_z 2 -stag_dof_2 2 -stag_grid_y 5
185 output_file: output/empty.out
186
187 test:
188 suffix: 3
189 nsize: 12
190 args: -stag_ranks_x 3 -stag_ranks_y 2 -stag_ranks_z 2 -stag_dof_0 2 -stag_grid_x 6
191 output_file: output/empty.out
192
193 test:
194 suffix: 4
195 nsize: 12
196 args: -stag_ranks_x 3 -stag_ranks_y 2 -stag_ranks_z 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_boundary_type_z ghosted -stag_stencil_width 1
197 output_file: output/empty.out
198
199 test:
200 suffix: 5
201 nsize: 12
202 args: -stag_ranks_x 3 -stag_ranks_y 2 -stag_ranks_z 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_z ghosted -stag_stencil_width 1
203 output_file: output/empty.out
204
205 test:
206 suffix: 6
207 nsize: 8
208 args: -stag_dof_0 3 -stag_dof_1 2 -stag_dof_2 4 -stag_dof_3 2 -stag_boundary_type_y ghosted -stag_boundary_type_z ghosted -stag_stencil_width 1
209 output_file: output/empty.out
210 TEST*/
211