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