xref: /petsc/src/snes/tests/ex2.c (revision 8fb5bd83c3955fefcf33a54e3bb66920a9fa884b) 
1 static char help[] = "Interpolation Tests for Plex\n\n";
2 
3 #include <petscsnes.h>
4 #include <petscdmplex.h>
5 #include <petscdmda.h>
6 #include <petscds.h>
7 
8 typedef enum {CENTROID, GRID, GRID_REPLICATED} PointType;
9 
10 typedef struct {
11   PointType pointType; /* Point generation mechanism */
12 } AppCtx;
13 
14 static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
15 {
16   PetscInt d, c;
17 
18   PetscCheck(Nc == 3,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Something is wrong: %" PetscInt_FMT, Nc);
19   for (c = 0; c < Nc; ++c) {
20     u[c] = 0.0;
21     for (d = 0; d < dim; ++d) u[c] += x[d];
22   }
23   return 0;
24 }
25 
26 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
27 {
28   const char    *pointTypes[3] = {"centroid", "grid", "grid_replicated"};
29   PetscInt       pt;
30 
31   PetscFunctionBegin;
32   options->pointType = CENTROID;
33   PetscOptionsBegin(comm, "", "Interpolation Options", "DMPLEX");
34   pt   = options->pointType;
35   PetscCall(PetscOptionsEList("-point_type", "The point type", "ex2.c", pointTypes, 3, pointTypes[options->pointType], &pt, NULL));
36   options->pointType = (PointType) pt;
37   PetscOptionsEnd();
38   PetscFunctionReturn(0);
39 }
40 
41 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm)
42 {
43   PetscFunctionBegin;
44   PetscCall(DMCreate(comm, dm));
45   PetscCall(DMSetType(*dm, DMPLEX));
46   PetscCall(DMSetFromOptions(*dm));
47   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
48   PetscFunctionReturn(0);
49 }
50 
51 static PetscErrorCode CreatePoints_Centroid(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx)
52 {
53   PetscSection   coordSection;
54   Vec            coordsLocal;
55   PetscInt       spaceDim, p;
56   PetscMPIInt    rank;
57 
58   PetscFunctionBegin;
59   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
60   PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
61   PetscCall(DMGetCoordinateSection(dm, &coordSection));
62   PetscCall(DMGetCoordinateDim(dm, &spaceDim));
63   PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, Np));
64   PetscCall(PetscCalloc1(*Np * spaceDim, pcoords));
65   for (p = 0; p < *Np; ++p) {
66     PetscScalar *coords = NULL;
67     PetscInt     size, num, n, d;
68 
69     PetscCall(DMPlexVecGetClosure(dm, coordSection, coordsLocal, p, &size, &coords));
70     num  = size/spaceDim;
71     for (n = 0; n < num; ++n) {
72       for (d = 0; d < spaceDim; ++d) (*pcoords)[p*spaceDim+d] += PetscRealPart(coords[n*spaceDim+d]) / num;
73     }
74     PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %" PetscInt_FMT " (", rank, p));
75     for (d = 0; d < spaceDim; ++d) {
76       PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%g", (double)(*pcoords)[p*spaceDim+d]));
77       if (d < spaceDim-1) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ", "));
78     }
79     PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ")\n"));
80     PetscCall(DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, p, &num, &coords));
81   }
82   PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL));
83   *pointsAllProcs = PETSC_FALSE;
84   PetscFunctionReturn(0);
85 }
86 
87 static PetscErrorCode CreatePoints_Grid(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx)
88 {
89   DM             da;
90   DMDALocalInfo  info;
91   PetscInt       N = 3, n = 0, dim, spaceDim, i, j, k, *ind, d;
92   PetscReal      *h;
93   PetscMPIInt    rank;
94 
95   PetscFunctionBegin;
96   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
97   PetscCall(DMGetDimension(dm, &dim));
98   PetscCall(DMGetCoordinateDim(dm, &spaceDim));
99   PetscCall(PetscCalloc1(spaceDim,&ind));
100   PetscCall(PetscCalloc1(spaceDim,&h));
101   h[0] = 1.0/(N-1); h[1] = 1.0/(N-1); h[2] = 1.0/(N-1);
102   PetscCall(DMDACreate(PetscObjectComm((PetscObject) dm), &da));
103   PetscCall(DMSetDimension(da, dim));
104   PetscCall(DMDASetSizes(da, N, N, N));
105   PetscCall(DMDASetDof(da, 1));
106   PetscCall(DMDASetStencilWidth(da, 1));
107   PetscCall(DMSetUp(da));
108   PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
109   PetscCall(DMDAGetLocalInfo(da, &info));
110   *Np  = info.xm * info.ym * info.zm;
111   PetscCall(PetscCalloc1(*Np * spaceDim, pcoords));
112   for (k = info.zs; k < info.zs + info.zm; ++k) {
113     ind[2] = k;
114     for (j = info.ys; j < info.ys + info.ym; ++j) {
115       ind[1] = j;
116       for (i = info.xs; i < info.xs + info.xm; ++i, ++n) {
117         ind[0] = i;
118 
119         for (d = 0; d < spaceDim; ++d) (*pcoords)[n*spaceDim+d] = ind[d]*h[d];
120         PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %" PetscInt_FMT " (", rank, n));
121         for (d = 0; d < spaceDim; ++d) {
122           PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%g", (double)(*pcoords)[n*spaceDim+d]));
123           if (d < spaceDim-1) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ", "));
124         }
125         PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ")\n"));
126       }
127     }
128   }
129   PetscCall(DMDestroy(&da));
130   PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL));
131   PetscCall(PetscFree(ind));
132   PetscCall(PetscFree(h));
133   *pointsAllProcs = PETSC_FALSE;
134   PetscFunctionReturn(0);
135 }
136 
137 static PetscErrorCode CreatePoints_GridReplicated(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx)
138 {
139   PetscInt       N = 3, n = 0, dim, spaceDim, i, j, k, *ind, d;
140   PetscReal      *h;
141   PetscMPIInt    rank;
142 
143   PetscFunctionBeginUser;
144   PetscCall(DMGetDimension(dm, &dim));
145   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
146   PetscCall(DMGetCoordinateDim(dm, &spaceDim));
147   PetscCall(PetscCalloc1(spaceDim,&ind));
148   PetscCall(PetscCalloc1(spaceDim,&h));
149   h[0] = 1.0/(N-1); h[1] = 1.0/(N-1); h[2] = 1.0/(N-1);
150   *Np  = N * (dim > 1 ? N : 1) * (dim > 2 ? N : 1);
151   PetscCall(PetscCalloc1(*Np * spaceDim, pcoords));
152   for (k = 0; k < N; ++k) {
153     ind[2] = k;
154     for (j = 0; j < N; ++j) {
155       ind[1] = j;
156       for (i = 0; i < N; ++i, ++n) {
157         ind[0] = i;
158 
159         for (d = 0; d < spaceDim; ++d) (*pcoords)[n*spaceDim+d] = ind[d]*h[d];
160         PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %" PetscInt_FMT " (", rank, n));
161         for (d = 0; d < spaceDim; ++d) {
162           PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%g", (double)(*pcoords)[n*spaceDim+d]));
163           if (d < spaceDim-1) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ", "));
164         }
165         PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ")\n"));
166       }
167     }
168   }
169   PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL));
170   *pointsAllProcs = PETSC_TRUE;
171   PetscCall(PetscFree(ind));
172   PetscCall(PetscFree(h));
173   PetscFunctionReturn(0);
174 }
175 
176 static PetscErrorCode CreatePoints(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx)
177 {
178   PetscFunctionBegin;
179   *pointsAllProcs = PETSC_FALSE;
180   switch (ctx->pointType) {
181   case CENTROID:        PetscCall(CreatePoints_Centroid(dm, Np, pcoords, pointsAllProcs, ctx));break;
182   case GRID:            PetscCall(CreatePoints_Grid(dm, Np, pcoords, pointsAllProcs, ctx));break;
183   case GRID_REPLICATED: PetscCall(CreatePoints_GridReplicated(dm, Np, pcoords, pointsAllProcs, ctx));break;
184   default: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Invalid point generation type %d", (int) ctx->pointType);
185   }
186   PetscFunctionReturn(0);
187 }
188 
189 int main(int argc, char **argv)
190 {
191   AppCtx              ctx;
192   PetscErrorCode   (**funcs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
193   DM                  dm;
194   PetscFE             fe;
195   DMInterpolationInfo interpolator;
196   Vec                 lu, fieldVals;
197   PetscScalar        *vals;
198   const PetscScalar  *ivals, *vcoords;
199   PetscReal          *pcoords;
200   PetscBool           simplex, pointsAllProcs=PETSC_TRUE;
201   PetscInt            dim, spaceDim, Nc, c, Np, p;
202   PetscMPIInt         rank, size;
203   PetscViewer         selfviewer;
204 
205   PetscCall(PetscInitialize(&argc, &argv, NULL,help));
206   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &ctx));
207   PetscCall(CreateMesh(PETSC_COMM_WORLD, &ctx, &dm));
208   PetscCall(DMGetDimension(dm, &dim));
209   PetscCall(DMGetCoordinateDim(dm, &spaceDim));
210   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
211   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
212   /* Create points */
213   PetscCall(CreatePoints(dm, &Np, &pcoords, &pointsAllProcs, &ctx));
214   /* Create interpolator */
215   PetscCall(DMInterpolationCreate(PETSC_COMM_WORLD, &interpolator));
216   PetscCall(DMInterpolationSetDim(interpolator, spaceDim));
217   PetscCall(DMInterpolationAddPoints(interpolator, Np, pcoords));
218   PetscCall(DMInterpolationSetUp(interpolator, dm, pointsAllProcs, PETSC_FALSE));
219   /* Check locations */
220   for (c = 0; c < interpolator->n; ++c) {
221     PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %" PetscInt_FMT " is in Cell %" PetscInt_FMT "\n", rank, c, interpolator->cells[c]));
222   }
223   PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL));
224   PetscCall(VecView(interpolator->coords, PETSC_VIEWER_STDOUT_WORLD));
225   /* Setup Discretization */
226   Nc   = dim;
227   PetscCall(DMPlexIsSimplex(dm, &simplex));
228   PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, Nc, simplex, NULL, -1, &fe));
229   PetscCall(DMSetField(dm, 0, NULL, (PetscObject) fe));
230   PetscCall(DMCreateDS(dm));
231   PetscCall(PetscFEDestroy(&fe));
232   /* Create function */
233   PetscCall(PetscCalloc2(Nc, &funcs, Nc, &vals));
234   for (c = 0; c < Nc; ++c) funcs[c] = linear;
235   PetscCall(DMGetLocalVector(dm, &lu));
236   PetscCall(DMProjectFunctionLocal(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, lu));
237   PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
238   PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD,PETSC_COMM_SELF,&selfviewer));
239   PetscCall(PetscViewerASCIIPrintf(selfviewer, "[%d]solution\n", rank));
240   PetscCall(VecView(lu,selfviewer));
241   PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD,PETSC_COMM_SELF,&selfviewer));
242   PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
243   PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));
244   /* Check interpolant */
245   PetscCall(VecCreateSeq(PETSC_COMM_SELF, interpolator->n * Nc, &fieldVals));
246   PetscCall(DMInterpolationSetDof(interpolator, Nc));
247   PetscCall(DMInterpolationEvaluate(interpolator, dm, lu, fieldVals));
248   for (p = 0; p < size; ++p) {
249     if (p == rank) {
250       PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d]Field values\n", rank));
251       PetscCall(VecView(fieldVals, PETSC_VIEWER_STDOUT_SELF));
252     }
253     PetscCall(PetscBarrier((PetscObject) dm));
254   }
255   PetscCall(VecGetArrayRead(interpolator->coords, &vcoords));
256   PetscCall(VecGetArrayRead(fieldVals, &ivals));
257   for (p = 0; p < interpolator->n; ++p) {
258     for (c = 0; c < Nc; ++c) {
259 #if defined(PETSC_USE_COMPLEX)
260       PetscReal vcoordsReal[3];
261       PetscInt  i;
262 
263       for (i = 0; i < spaceDim; i++) vcoordsReal[i] = PetscRealPart(vcoords[p * spaceDim + i]);
264 #else
265       const PetscReal *vcoordsReal = &vcoords[p*spaceDim];
266 #endif
267       (*funcs[c])(dim, 0.0, vcoordsReal, Nc, vals, NULL);
268       if (PetscAbsScalar(ivals[p*Nc+c] - vals[c]) > PETSC_SQRT_MACHINE_EPSILON)
269         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid interpolated value %g != %g (%" PetscInt_FMT ", %" PetscInt_FMT ")", (double) PetscRealPart(ivals[p*Nc+c]), (double) PetscRealPart(vals[c]), p, c);
270     }
271   }
272   PetscCall(VecRestoreArrayRead(interpolator->coords, &vcoords));
273   PetscCall(VecRestoreArrayRead(fieldVals, &ivals));
274   /* Cleanup */
275   PetscCall(PetscFree(pcoords));
276   PetscCall(PetscFree2(funcs, vals));
277   PetscCall(VecDestroy(&fieldVals));
278   PetscCall(DMRestoreLocalVector(dm, &lu));
279   PetscCall(DMInterpolationDestroy(&interpolator));
280   PetscCall(DMDestroy(&dm));
281   PetscCall(PetscFinalize());
282   return 0;
283 }
284 
285 /*TEST
286 
287   testset:
288     requires: ctetgen
289     args: -dm_plex_dim 3 -petscspace_degree 1
290 
291     test:
292       suffix: 0
293     test:
294       suffix: 1
295       args: -dm_refine 2
296     test:
297       suffix: 2
298       nsize: 2
299       args: -petscpartitioner_type simple
300     test:
301       suffix: 3
302       nsize: 2
303       args: -dm_refine 2 -petscpartitioner_type simple
304     test:
305       suffix: 4
306       nsize: 5
307       args: -petscpartitioner_type simple
308     test:
309       suffix: 5
310       nsize: 5
311       args: -dm_refine 2 -petscpartitioner_type simple
312     test:
313       suffix: 6
314       args: -point_type grid
315     test:
316       suffix: 7
317       args: -dm_refine 2 -point_type grid
318     test:
319       suffix: 8
320       nsize: 2
321       args: -petscpartitioner_type simple -point_type grid
322     test:
323       suffix: 9
324       args: -point_type grid_replicated
325     test:
326       suffix: 10
327       nsize: 2
328       args: -petscpartitioner_type simple -point_type grid_replicated
329     test:
330       suffix: 11
331       nsize: 2
332       args: -dm_refine 2 -petscpartitioner_type simple -point_type grid_replicated
333 
334 TEST*/
335