1 static char help[] = "(Partially) test DMStag default interpolation, 2d faces-only.\n\n";
2
3 #include <petscdm.h>
4 #include <petscdmstag.h>
5 #include <petscksp.h>
6
7 PetscErrorCode CreateSystem(DM dm, Mat *A, Vec *b);
8
main(int argc,char ** argv)9 int main(int argc, char **argv)
10 {
11 DM dm, dmCoarse;
12 Mat Ai;
13
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));
21
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 }
78
79 PetscCall(MatDestroy(&Ai));
80 PetscCall(DMDestroy(&dm));
81 PetscCall(DMDestroy(&dmCoarse));
82 PetscCall(PetscFinalize());
83 return 0;
84 }
85
86 /*TEST
87
88 test:
89 suffix: 1
90 nsize: 1
91 args:
92 output_file: output/empty.out
93
94 test:
95 suffix: 2
96 nsize: 4
97 args: -stag_grid_x 8 -stag_grid_y 4
98 output_file: output/empty.out
99
100 TEST*/
101