Actual source code: ex7.c

  1: static char help[] = "Test DMStag 3d periodic and ghosted boundary conditions\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;
 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;

 15:   PetscInitialize(&argc, &argv, (char *)0, help);
 16:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 17:   dof0         = 1;
 18:   dof1         = 1;
 19:   dof2         = 1;
 20:   dof3         = 1;
 21:   stencilWidth = 2;
 22:   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:   DMSetFromOptions(dm);
 24:   DMSetUp(dm);
 25:   DMStagGetDOF(dm, &dof0, &dof1, &dof2, &dof3);
 26:   dofTotal = dof0 + 3 * dof1 + 3 * dof2 + dof3;
 27:   DMStagGetStencilWidth(dm, &stencilWidth);

 29:   DMCreateLocalVector(dm, &vecLocal1);
 30:   VecDuplicate(vecLocal1, &vecLocal2);

 32:   DMCreateGlobalVector(dm, &vec);
 33:   VecSet(vec, 1.0);
 34:   VecSet(vecLocal1, 0.0);
 35:   DMGlobalToLocalBegin(dm, vec, INSERT_VALUES, vecLocal1);
 36:   DMGlobalToLocalEnd(dm, vec, INSERT_VALUES, vecLocal1);

 38:   DMStagGetCorners(dm, &startx, &starty, &startz, &nx, &ny, &nz, NULL, NULL, NULL);
 39:   DMStagVecGetArrayRead(dm, vecLocal1, &a1);
 40:   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) 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:   DMStagVecRestoreArrayRead(dm, vecLocal1, &a1);
 57:   DMStagVecRestoreArray(dm, vecLocal2, &a2);

 59:   DMLocalToGlobalBegin(dm, vecLocal2, INSERT_VALUES, vec);
 60:   DMLocalToGlobalEnd(dm, vecLocal2, INSERT_VALUES, vec);

 62:   /* For the all-periodic case, all values are the same . Otherwise, just check the local version */
 63:   DMStagGetBoundaryTypes(dm, &boundaryTypex, &boundaryTypey, &boundaryTypez);
 64:   if (boundaryTypex == DM_BOUNDARY_PERIODIC && boundaryTypey == DM_BOUNDARY_PERIODIC && boundaryTypez == DM_BOUNDARY_PERIODIC) {
 65:     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) PetscPrintf(PETSC_COMM_SELF, "[%d] Unexpected value %g (expecting %g)\n", rank, (double)PetscRealPart(a[i]), (double)PetscRealPart(expected));
 70:     }
 71:     VecRestoreArray(vec, &a);
 72:   } else {
 73:     DMStagVecGetArrayRead(dm, vecLocal2, &a2);
 74:     DMStagGetGlobalSizes(dm, &Nx, &Ny, &Nz);
 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:                 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:                 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:                 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:                 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:                 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:                 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:                 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:                 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:     DMStagVecRestoreArrayRead(dm, vecLocal2, &a2);
163:   }

165:   VecDestroy(&vec);
166:   VecDestroy(&vecLocal1);
167:   VecDestroy(&vecLocal2);
168:   DMDestroy(&dm);
169:   PetscFinalize();
170:   return 0;
171: }

173: /*TEST

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

180:    test:
181:       suffix: 2
182:       nsize: 8
183:       args: -stag_ranks_x 2 -stag_ranks_y 2 -stag_ranks_z 2 -stag_dof_2 2 -stag_grid_y 5

185:    test:
186:       suffix: 3
187:       nsize: 12
188:       args: -stag_ranks_x 3 -stag_ranks_y 2 -stag_ranks_z 2 -stag_dof_0 2 -stag_grid_x 6

190:    test:
191:       suffix: 4
192:       nsize: 12
193:       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

195:    test:
196:       suffix: 5
197:       nsize: 12
198:       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

200:    test:
201:       suffix: 6
202:       nsize: 8
203:       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
204: TEST*/