Actual source code: ex10.c

  1: static char help[] = "Test DMStag 2d 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, nx, ny, i, j, d, is, js, dof0, dof1, dof2, dofTotal, stencilWidth, Nx, Ny;
 11:   DMBoundaryType boundaryTypex, boundaryTypey;
 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:   stencilWidth = 2;
 21:   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:   DMSetFromOptions(dm);
 23:   DMSetUp(dm);
 24:   DMStagGetDOF(dm, &dof0, &dof1, &dof2, NULL);
 25:   dofTotal = dof0 + 2 * dof1 + dof2;
 26:   DMStagGetStencilWidth(dm, &stencilWidth);

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

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

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

 54:   DMLocalToGlobalBegin(dm, vecLocal2, INSERT_VALUES, vec);
 55:   DMLocalToGlobalEnd(dm, vecLocal2, INSERT_VALUES, vec);

 57:   /* For the all-periodic case, all values are the same . Otherwise, just check the local version */
 58:   DMStagGetBoundaryTypes(dm, &boundaryTypex, &boundaryTypey, NULL);
 59:   if (boundaryTypex == DM_BOUNDARY_PERIODIC && boundaryTypey == DM_BOUNDARY_PERIODIC) {
 60:     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) PetscPrintf(PETSC_COMM_SELF, "[%d] Unexpected value %g (expecting %g)\n", rank, (double)PetscRealPart(a[i]), (double)PetscRealPart(expected));
 65:     }
 66:     VecRestoreArray(vec, &a);
 67:   } else {
 68:     DMStagVecGetArrayRead(dm, vecLocal2, &a2);
 69:     DMStagGetGlobalSizes(dm, &Nx, &Ny, NULL);
 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:               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:               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:               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:               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:     DMStagVecRestoreArrayRead(dm, vecLocal2, &a2);
118:   }

120:   VecDestroy(&vec);
121:   VecDestroy(&vecLocal1);
122:   VecDestroy(&vecLocal2);
123:   DMDestroy(&dm);
124:   PetscFinalize();
125:   return 0;
126: }

128: /*TEST

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

135:    test:
136:       suffix: 2
137:       nsize: 4
138:       args: -stag_ranks_x 2 -stag_ranks_y 2 -stag_dof_1 2 -stag_grid_y 5

140:    test:
141:       suffix: 3
142:       nsize: 6
143:       args: -stag_ranks_x 3 -stag_ranks_y 2 -stag_dof_0 2 -stag_grid_x 6

145:    test:
146:       suffix: 4
147:       nsize: 6
148:       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

150:    test:
151:       suffix: 5
152:       nsize: 6
153:       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

155:    test:
156:       suffix: 6
157:       nsize: 9
158:       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

160:    test:
161:       suffix: 7
162:       nsize: 1
163:       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

165: TEST*/