Actual source code: ex19.c

  1: static char help[] = "(Partially) test DMStag default interpolation, 2d faces-only.\n\n";

  3: #include <petscdm.h>
  4: #include <petscdmstag.h>
  5: #include <petscksp.h>

  7: PetscErrorCode CreateSystem(DM dm, Mat *A, Vec *b);

  9: int main(int argc, char **argv)
 10: {
 11:   DM  dm, dmCoarse;
 12:   Mat Ai;

 14:   PetscFunctionBeginUser;
 15:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 16:   PetscCall(DMStagCreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, 2, 4, PETSC_DECIDE, PETSC_DECIDE, 0, 1, 1, DMSTAG_STENCIL_BOX, 1, NULL, NULL, &dm));
 17:   PetscCall(DMSetFromOptions(dm));
 18:   PetscCall(DMSetUp(dm));
 19:   PetscCall(DMCoarsen(dm, MPI_COMM_NULL, &dmCoarse));
 20:   PetscCall(DMCreateInterpolation(dmCoarse, dm, &Ai, NULL));

 22:   /* See what happens to a constant value on each sub-grid */
 23:   {
 24:     Vec localCoarse, globalCoarse, globalFine, localFine;
 25:     PetscCall(DMGetGlobalVector(dm, &globalFine));
 26:     PetscCall(DMGetGlobalVector(dmCoarse, &globalCoarse));
 27:     PetscCall(DMGetLocalVector(dmCoarse, &localCoarse));
 28:     PetscCall(DMGetLocalVector(dm, &localFine));
 29:     PetscCall(VecSet(localCoarse, -1.0));
 30:     PetscCall(VecSet(localFine, -1.0));
 31:     {
 32:       PetscInt       i, j, startx, starty, nx, ny, extrax, extray;
 33:       PetscInt       p, vx, vy;
 34:       PetscScalar ***arr;
 35:       PetscCall(DMStagGetCorners(dmCoarse, &startx, &starty, NULL, &nx, &ny, NULL, &extrax, &extray, NULL));
 36:       PetscCall(DMStagVecGetArray(dmCoarse, localCoarse, &arr));
 37:       PetscCall(DMStagGetLocationSlot(dmCoarse, DMSTAG_LEFT, 0, &vx));
 38:       PetscCall(DMStagGetLocationSlot(dmCoarse, DMSTAG_DOWN, 0, &vy));
 39:       PetscCall(DMStagGetLocationSlot(dmCoarse, DMSTAG_ELEMENT, 0, &p));
 40:       for (j = starty; j < starty + ny + extray; ++j) {
 41:         for (i = startx; i < startx + nx + extrax; ++i) {
 42:           arr[j][i][vy] = (i < startx + nx) ? 10.0 : -1;
 43:           arr[j][i][vx] = (j < starty + ny) ? 20.0 : -1;
 44:           arr[j][i][p]  = (i < startx + nx) && (j < starty + ny) ? 30.0 : -1;
 45:         }
 46:       }
 47:       PetscCall(DMStagVecRestoreArray(dmCoarse, localCoarse, &arr));
 48:     }
 49:     PetscCall(DMLocalToGlobal(dmCoarse, localCoarse, INSERT_VALUES, globalCoarse));
 50:     PetscCall(MatInterpolate(Ai, globalCoarse, globalFine));
 51:     PetscCall(DMGlobalToLocal(dm, globalFine, INSERT_VALUES, localFine));
 52:     {
 53:       PetscInt       i, j, startx, starty, nx, ny, extrax, extray;
 54:       PetscInt       p, vx, vy;
 55:       PetscScalar ***arr;
 56:       PetscCall(DMStagGetCorners(dm, &startx, &starty, NULL, &nx, &ny, NULL, &extrax, &extray, NULL));
 57:       PetscCall(DMStagVecGetArrayRead(dm, localFine, &arr));
 58:       PetscCall(DMStagGetLocationSlot(dm, DMSTAG_LEFT, 0, &vx));
 59:       PetscCall(DMStagGetLocationSlot(dm, DMSTAG_DOWN, 0, &vy));
 60:       PetscCall(DMStagGetLocationSlot(dm, DMSTAG_ELEMENT, 0, &p));
 61:       for (j = starty; j < starty + ny + extray; ++j) {
 62:         for (i = startx; i < startx + nx + extrax; ++i) {
 63:           const PetscScalar expected_vy = (i < startx + nx) ? 10.0 : -1;
 64:           const PetscScalar expected_vx = (j < starty + ny) ? 20.0 : -1;
 65:           const PetscScalar expected_p  = (i < startx + nx) && (j < starty + ny) ? 30.0 : -1;
 66:           if (arr[j][i][vy] != expected_vy) PetscCall(PetscPrintf(PETSC_COMM_SELF, "wrong %" PetscInt_FMT " %" PetscInt_FMT "\n", i, j));
 67:           if (arr[j][i][vx] != expected_vx) PetscCall(PetscPrintf(PETSC_COMM_SELF, "wrong %" PetscInt_FMT " %" PetscInt_FMT "\n", i, j));
 68:           if (arr[j][i][p] != expected_p) PetscCall(PetscPrintf(PETSC_COMM_SELF, "wrong %" PetscInt_FMT " %" PetscInt_FMT "\n", i, j));
 69:         }
 70:       }
 71:       PetscCall(DMStagVecRestoreArrayRead(dm, localFine, &arr));
 72:     }
 73:     PetscCall(DMRestoreLocalVector(dmCoarse, &localCoarse));
 74:     PetscCall(DMRestoreLocalVector(dm, &localFine));
 75:     PetscCall(DMRestoreGlobalVector(dmCoarse, &globalCoarse));
 76:     PetscCall(DMRestoreGlobalVector(dm, &globalFine));
 77:   }

 79:   PetscCall(MatDestroy(&Ai));
 80:   PetscCall(DMDestroy(&dm));
 81:   PetscCall(DMDestroy(&dmCoarse));
 82:   PetscCall(PetscFinalize());
 83:   return 0;
 84: }

 86: /*TEST

 88:    test:
 89:       suffix: 1
 90:       nsize: 1
 91:       args:

 93:    test:
 94:       suffix: 2
 95:       nsize: 4
 96:       args: -stag_grid_x 8 -stag_grid_y 4

 98: TEST*/