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;

 15:   PetscInitialize(&argc, &argv, (char *)0, help);
 16:   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:   DMSetFromOptions(dm);
 18:   DMSetUp(dm);
 19:   DMCoarsen(dm, MPI_COMM_NULL, &dmCoarse);
 20:   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:     DMGetGlobalVector(dm, &globalFine);
 26:     DMGetGlobalVector(dmCoarse, &globalCoarse);
 27:     DMGetLocalVector(dmCoarse, &localCoarse);
 28:     DMGetLocalVector(dm, &localFine);
 29:     VecSet(localCoarse, -1.0);
 30:     VecSet(localFine, -1.0);
 31:     {
 32:       PetscInt       i, j, startx, starty, nx, ny, extrax, extray;
 33:       PetscInt       p, vx, vy;
 34:       PetscScalar ***arr;
 35:       DMStagGetCorners(dmCoarse, &startx, &starty, NULL, &nx, &ny, NULL, &extrax, &extray, NULL);
 36:       DMStagVecGetArray(dmCoarse, localCoarse, &arr);
 37:       DMStagGetLocationSlot(dmCoarse, DMSTAG_LEFT, 0, &vx);
 38:       DMStagGetLocationSlot(dmCoarse, DMSTAG_DOWN, 0, &vy);
 39:       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:       DMStagVecRestoreArray(dmCoarse, localCoarse, &arr);
 48:     }
 49:     DMLocalToGlobal(dmCoarse, localCoarse, INSERT_VALUES, globalCoarse);
 50:     MatInterpolate(Ai, globalCoarse, globalFine);
 51:     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:       DMStagGetCorners(dm, &startx, &starty, NULL, &nx, &ny, NULL, &extrax, &extray, NULL);
 57:       DMStagVecGetArrayRead(dm, localFine, &arr);
 58:       DMStagGetLocationSlot(dm, DMSTAG_LEFT, 0, &vx);
 59:       DMStagGetLocationSlot(dm, DMSTAG_DOWN, 0, &vy);
 60:       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) PetscPrintf(PETSC_COMM_SELF, "wrong %" PetscInt_FMT " %" PetscInt_FMT "\n", i, j);
 67:           if (arr[j][i][vx] != expected_vx) PetscPrintf(PETSC_COMM_SELF, "wrong %" PetscInt_FMT " %" PetscInt_FMT "\n", i, j);
 68:           if (arr[j][i][p] != expected_p) PetscPrintf(PETSC_COMM_SELF, "wrong %" PetscInt_FMT " %" PetscInt_FMT "\n", i, j);
 69:         }
 70:       }
 71:       DMStagVecRestoreArrayRead(dm, localFine, &arr);
 72:     }
 73:     DMRestoreLocalVector(dmCoarse, &localCoarse);
 74:     DMRestoreLocalVector(dm, &localFine);
 75:     DMRestoreGlobalVector(dmCoarse, &globalCoarse);
 76:     DMRestoreGlobalVector(dm, &globalFine);
 77:   }

 79:   MatDestroy(&Ai);
 80:   DMDestroy(&dm);
 81:   DMDestroy(&dmCoarse);
 82:   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*/