xref: /petsc/src/snes/tests/ex2.c (revision 800f99ff9e85495c69e9e5819c0be0dbd8cbc57c) 
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   PetscFunctionBeginUser;
206   PetscCall(PetscInitialize(&argc, &argv, NULL,help));
207   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &ctx));
208   PetscCall(CreateMesh(PETSC_COMM_WORLD, &ctx, &dm));
209   PetscCall(DMGetDimension(dm, &dim));
210   PetscCall(DMGetCoordinateDim(dm, &spaceDim));
211   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
212   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
213   /* Create points */
214   PetscCall(CreatePoints(dm, &Np, &pcoords, &pointsAllProcs, &ctx));
215   /* Create interpolator */
216   PetscCall(DMInterpolationCreate(PETSC_COMM_WORLD, &interpolator));
217   PetscCall(DMInterpolationSetDim(interpolator, spaceDim));
218   PetscCall(DMInterpolationAddPoints(interpolator, Np, pcoords));
219   PetscCall(DMInterpolationSetUp(interpolator, dm, pointsAllProcs, PETSC_FALSE));
220   /* Check locations */
221   for (c = 0; c < interpolator->n; ++c) {
222     PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %" PetscInt_FMT " is in Cell %" PetscInt_FMT "\n", rank, c, interpolator->cells[c]));
223   }
224   PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL));
225   PetscCall(VecView(interpolator->coords, PETSC_VIEWER_STDOUT_WORLD));
226   /* Setup Discretization */
227   Nc   = dim;
228   PetscCall(DMPlexIsSimplex(dm, &simplex));
229   PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, Nc, simplex, NULL, -1, &fe));
230   PetscCall(DMSetField(dm, 0, NULL, (PetscObject) fe));
231   PetscCall(DMCreateDS(dm));
232   PetscCall(PetscFEDestroy(&fe));
233   /* Create function */
234   PetscCall(PetscCalloc2(Nc, &funcs, Nc, &vals));
235   for (c = 0; c < Nc; ++c) funcs[c] = linear;
236   PetscCall(DMGetLocalVector(dm, &lu));
237   PetscCall(DMProjectFunctionLocal(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, lu));
238   PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
239   PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD,PETSC_COMM_SELF,&selfviewer));
240   PetscCall(PetscViewerASCIIPrintf(selfviewer, "[%d]solution\n", rank));
241   PetscCall(VecView(lu,selfviewer));
242   PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD,PETSC_COMM_SELF,&selfviewer));
243   PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
244   PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));
245   /* Check interpolant */
246   PetscCall(VecCreateSeq(PETSC_COMM_SELF, interpolator->n * Nc, &fieldVals));
247   PetscCall(DMInterpolationSetDof(interpolator, Nc));
248   PetscCall(DMInterpolationEvaluate(interpolator, dm, lu, fieldVals));
249   for (p = 0; p < size; ++p) {
250     if (p == rank) {
251       PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d]Field values\n", rank));
252       PetscCall(VecView(fieldVals, PETSC_VIEWER_STDOUT_SELF));
253     }
254     PetscCall(PetscBarrier((PetscObject) dm));
255   }
256   PetscCall(VecGetArrayRead(interpolator->coords, &vcoords));
257   PetscCall(VecGetArrayRead(fieldVals, &ivals));
258   for (p = 0; p < interpolator->n; ++p) {
259     for (c = 0; c < Nc; ++c) {
260 #if defined(PETSC_USE_COMPLEX)
261       PetscReal vcoordsReal[3];
262       PetscInt  i;
263 
264       for (i = 0; i < spaceDim; i++) vcoordsReal[i] = PetscRealPart(vcoords[p * spaceDim + i]);
265 #else
266       const PetscReal *vcoordsReal = &vcoords[p*spaceDim];
267 #endif
268       (*funcs[c])(dim, 0.0, vcoordsReal, Nc, vals, NULL);
269       if (PetscAbsScalar(ivals[p*Nc+c] - vals[c]) > PETSC_SQRT_MACHINE_EPSILON)
270         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);
271     }
272   }
273   PetscCall(VecRestoreArrayRead(interpolator->coords, &vcoords));
274   PetscCall(VecRestoreArrayRead(fieldVals, &ivals));
275   /* Cleanup */
276   PetscCall(PetscFree(pcoords));
277   PetscCall(PetscFree2(funcs, vals));
278   PetscCall(VecDestroy(&fieldVals));
279   PetscCall(DMRestoreLocalVector(dm, &lu));
280   PetscCall(DMInterpolationDestroy(&interpolator));
281   PetscCall(DMDestroy(&dm));
282   PetscCall(PetscFinalize());
283   return 0;
284 }
285 
286 /*TEST
287 
288   testset:
289     requires: ctetgen
290     args: -dm_plex_dim 3 -petscspace_degree 1
291 
292     test:
293       suffix: 0
294     test:
295       suffix: 1
296       args: -dm_refine 2
297     test:
298       suffix: 2
299       nsize: 2
300       args: -petscpartitioner_type simple
301     test:
302       suffix: 3
303       nsize: 2
304       args: -dm_refine 2 -petscpartitioner_type simple
305     test:
306       suffix: 4
307       nsize: 5
308       args: -petscpartitioner_type simple
309     test:
310       suffix: 5
311       nsize: 5
312       args: -dm_refine 2 -petscpartitioner_type simple
313     test:
314       suffix: 6
315       args: -point_type grid
316     test:
317       suffix: 7
318       args: -dm_refine 2 -point_type grid
319     test:
320       suffix: 8
321       nsize: 2
322       args: -petscpartitioner_type simple -point_type grid
323     test:
324       suffix: 9
325       args: -point_type grid_replicated
326     test:
327       suffix: 10
328       nsize: 2
329       args: -petscpartitioner_type simple -point_type grid_replicated
330     test:
331       suffix: 11
332       nsize: 2
333       args: -dm_refine 2 -petscpartitioner_type simple -point_type grid_replicated
334 
335 TEST*/
336