Actual source code: ex2.c

  1: static char help[] = "Interpolation Tests for Plex\n\n";

  3: #include <petscsnes.h>
  4: #include <petscdmplex.h>
  5: #include <petscdmda.h>
  6: #include <petscds.h>

  8: typedef enum {
  9:   CENTROID,
 10:   GRID,
 11:   GRID_REPLICATED
 12: } PointType;

 14: typedef enum {
 15:   CONSTANT,
 16:   LINEAR
 17: } FuncType;

 19: typedef struct {
 20:   PointType pointType; // Point generation mechanism
 21:   FuncType  funcType;  // Type of interpolated function
 22:   PetscBool useFV;     // Use finite volume, instead of finite element
 23: } AppCtx;

 25: static PetscErrorCode constant(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 26: {
 27:   PetscFunctionBeginUser;
 28:   for (PetscInt c = 0; c < Nc; ++c) u[c] = c + 1.;
 29:   PetscFunctionReturn(PETSC_SUCCESS);
 30: }

 32: static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 33: {
 34:   PetscFunctionBeginUser;
 35:   PetscCheck(Nc == 3, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Something is wrong: %" PetscInt_FMT, Nc);
 36:   for (PetscInt c = 0; c < Nc; ++c) {
 37:     u[c] = 0.0;
 38:     for (PetscInt d = 0; d < dim; ++d) u[c] += x[d];
 39:   }
 40:   PetscFunctionReturn(PETSC_SUCCESS);
 41: }

 43: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 44: {
 45:   const char *pointTypes[3] = {"centroid", "grid", "grid_replicated"};
 46:   const char *funcTypes[2]  = {"constant", "linear"};
 47:   PetscInt    pt, fn;

 49:   PetscFunctionBegin;
 50:   options->pointType = CENTROID;
 51:   options->funcType  = LINEAR;
 52:   options->useFV     = PETSC_FALSE;

 54:   PetscOptionsBegin(comm, "", "Interpolation Options", "DMPLEX");
 55:   pt = options->pointType;
 56:   PetscCall(PetscOptionsEList("-point_type", "The point type", "ex2.c", pointTypes, 3, pointTypes[options->pointType], &pt, NULL));
 57:   options->pointType = (PointType)pt;
 58:   fn                 = options->funcType;
 59:   PetscCall(PetscOptionsEList("-func_type", "The function type", "ex2.c", funcTypes, 2, funcTypes[options->funcType], &fn, NULL));
 60:   options->funcType = (FuncType)fn;
 61:   PetscCall(PetscOptionsBool("-use_fv", "Use finite volumes, instead of finite elements", "ex2.c", options->useFV, &options->useFV, NULL));
 62:   PetscOptionsEnd();
 63:   PetscFunctionReturn(PETSC_SUCCESS);
 64: }

 66: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm)
 67: {
 68:   PetscFunctionBegin;
 69:   PetscCall(DMCreate(comm, dm));
 70:   PetscCall(DMSetType(*dm, DMPLEX));
 71:   PetscCall(DMSetFromOptions(*dm));
 72:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
 73:   PetscFunctionReturn(PETSC_SUCCESS);
 74: }

 76: static PetscErrorCode CreatePoints_Centroid(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx)
 77: {
 78:   PetscSection coordSection;
 79:   Vec          coordsLocal;
 80:   PetscInt     spaceDim, p;
 81:   PetscMPIInt  rank;

 83:   PetscFunctionBegin;
 84:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 85:   PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
 86:   PetscCall(DMGetCoordinateSection(dm, &coordSection));
 87:   PetscCall(DMGetCoordinateDim(dm, &spaceDim));
 88:   PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, Np));
 89:   PetscCall(PetscCalloc1(*Np * spaceDim, pcoords));
 90:   for (p = 0; p < *Np; ++p) {
 91:     PetscScalar *coords = NULL;
 92:     PetscInt     size, num, n, d;

 94:     PetscCall(DMPlexVecGetClosure(dm, coordSection, coordsLocal, p, &size, &coords));
 95:     num = size / spaceDim;
 96:     for (n = 0; n < num; ++n) {
 97:       for (d = 0; d < spaceDim; ++d) (*pcoords)[p * spaceDim + d] += PetscRealPart(coords[n * spaceDim + d]) / num;
 98:     }
 99:     PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %" PetscInt_FMT " (", rank, p));
100:     for (d = 0; d < spaceDim; ++d) {
101:       PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%g", (double)(*pcoords)[p * spaceDim + d]));
102:       if (d < spaceDim - 1) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ", "));
103:     }
104:     PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ")\n"));
105:     PetscCall(DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, p, &num, &coords));
106:   }
107:   PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL));
108:   *pointsAllProcs = PETSC_FALSE;
109:   PetscFunctionReturn(PETSC_SUCCESS);
110: }

112: static PetscErrorCode CreatePoints_Grid(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx)
113: {
114:   DM            da;
115:   DMDALocalInfo info;
116:   PetscInt      N = 3, n = 0, dim, spaceDim, i, j, k, *ind, d;
117:   PetscReal    *h;
118:   PetscMPIInt   rank;

120:   PetscFunctionBegin;
121:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
122:   PetscCall(DMGetDimension(dm, &dim));
123:   PetscCall(DMGetCoordinateDim(dm, &spaceDim));
124:   PetscCall(PetscCalloc1(spaceDim, &ind));
125:   PetscCall(PetscCalloc1(spaceDim, &h));
126:   h[0] = 1.0 / (N - 1);
127:   h[1] = 1.0 / (N - 1);
128:   h[2] = 1.0 / (N - 1);
129:   PetscCall(DMDACreate(PetscObjectComm((PetscObject)dm), &da));
130:   PetscCall(DMSetDimension(da, dim));
131:   PetscCall(DMDASetSizes(da, N, N, N));
132:   PetscCall(DMDASetDof(da, 1));
133:   PetscCall(DMDASetStencilWidth(da, 1));
134:   PetscCall(DMSetUp(da));
135:   PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
136:   PetscCall(DMDAGetLocalInfo(da, &info));
137:   *Np = info.xm * info.ym * info.zm;
138:   PetscCall(PetscCalloc1(*Np * spaceDim, pcoords));
139:   for (k = info.zs; k < info.zs + info.zm; ++k) {
140:     ind[2] = k;
141:     for (j = info.ys; j < info.ys + info.ym; ++j) {
142:       ind[1] = j;
143:       for (i = info.xs; i < info.xs + info.xm; ++i, ++n) {
144:         ind[0] = i;

146:         for (d = 0; d < spaceDim; ++d) (*pcoords)[n * spaceDim + d] = ind[d] * h[d];
147:         PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %" PetscInt_FMT " (", rank, n));
148:         for (d = 0; d < spaceDim; ++d) {
149:           PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%g", (double)(*pcoords)[n * spaceDim + d]));
150:           if (d < spaceDim - 1) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ", "));
151:         }
152:         PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ")\n"));
153:       }
154:     }
155:   }
156:   PetscCall(DMDestroy(&da));
157:   PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL));
158:   PetscCall(PetscFree(ind));
159:   PetscCall(PetscFree(h));
160:   *pointsAllProcs = PETSC_FALSE;
161:   PetscFunctionReturn(PETSC_SUCCESS);
162: }

164: static PetscErrorCode CreatePoints_GridReplicated(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx)
165: {
166:   PetscInt    N = 3, n = 0, dim, spaceDim, i, j, k, *ind, d;
167:   PetscReal  *h;
168:   PetscMPIInt rank;

170:   PetscFunctionBeginUser;
171:   PetscCall(DMGetDimension(dm, &dim));
172:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
173:   PetscCall(DMGetCoordinateDim(dm, &spaceDim));
174:   PetscCall(PetscCalloc1(spaceDim, &ind));
175:   PetscCall(PetscCalloc1(spaceDim, &h));
176:   h[0] = 1.0 / (N - 1);
177:   h[1] = 1.0 / (N - 1);
178:   h[2] = 1.0 / (N - 1);
179:   *Np  = N * (dim > 1 ? N : 1) * (dim > 2 ? N : 1);
180:   PetscCall(PetscCalloc1(*Np * spaceDim, pcoords));
181:   for (k = 0; k < N; ++k) {
182:     ind[2] = k;
183:     for (j = 0; j < N; ++j) {
184:       ind[1] = j;
185:       for (i = 0; i < N; ++i, ++n) {
186:         ind[0] = i;

188:         for (d = 0; d < spaceDim; ++d) (*pcoords)[n * spaceDim + d] = ind[d] * h[d];
189:         PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %" PetscInt_FMT " (", rank, n));
190:         for (d = 0; d < spaceDim; ++d) {
191:           PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%g", (double)(*pcoords)[n * spaceDim + d]));
192:           if (d < spaceDim - 1) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ", "));
193:         }
194:         PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ")\n"));
195:       }
196:     }
197:   }
198:   PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL));
199:   *pointsAllProcs = PETSC_TRUE;
200:   PetscCall(PetscFree(ind));
201:   PetscCall(PetscFree(h));
202:   PetscFunctionReturn(PETSC_SUCCESS);
203: }

205: static PetscErrorCode CreatePoints(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx)
206: {
207:   PetscFunctionBegin;
208:   *pointsAllProcs = PETSC_FALSE;
209:   switch (ctx->pointType) {
210:   case CENTROID:
211:     PetscCall(CreatePoints_Centroid(dm, Np, pcoords, pointsAllProcs, ctx));
212:     break;
213:   case GRID:
214:     PetscCall(CreatePoints_Grid(dm, Np, pcoords, pointsAllProcs, ctx));
215:     break;
216:   case GRID_REPLICATED:
217:     PetscCall(CreatePoints_GridReplicated(dm, Np, pcoords, pointsAllProcs, ctx));
218:     break;
219:   default:
220:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid point generation type %d", (int)ctx->pointType);
221:   }
222:   PetscFunctionReturn(PETSC_SUCCESS);
223: }

225: static PetscErrorCode CreateDiscretization(DM dm, PetscInt Nc, AppCtx *ctx)
226: {
227:   PetscFunctionBegin;
228:   if (ctx->useFV) {
229:     PetscFV  fv;
230:     PetscInt cdim;

232:     PetscCall(PetscFVCreate(PetscObjectComm((PetscObject)dm), &fv));
233:     PetscCall(PetscObjectSetName((PetscObject)fv, "phi"));
234:     PetscCall(PetscFVSetFromOptions(fv));
235:     PetscCall(PetscFVSetNumComponents(fv, Nc));
236:     PetscCall(DMGetCoordinateDim(dm, &cdim));
237:     PetscCall(PetscFVSetSpatialDimension(fv, cdim));
238:     PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fv));
239:     PetscCall(PetscFVDestroy(&fv));
240:     PetscCall(DMCreateDS(dm));
241:   } else {
242:     PetscFE        fe;
243:     DMPolytopeType ct;
244:     PetscInt       dim, cStart;

246:     PetscCall(DMGetDimension(dm, &dim));
247:     PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
248:     PetscCall(DMPlexGetCellType(dm, cStart, &ct));
249:     PetscCall(PetscFECreateByCell(PetscObjectComm((PetscObject)dm), dim, Nc, ct, NULL, -1, &fe));
250:     PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
251:     PetscCall(PetscFEDestroy(&fe));
252:     PetscCall(DMCreateDS(dm));
253:   }
254:   PetscFunctionReturn(PETSC_SUCCESS);
255: }

257: int main(int argc, char **argv)
258: {
259:   AppCtx ctx;
260:   PetscErrorCode (**funcs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
261:   DM                  dm;
262:   DMInterpolationInfo interpolator;
263:   Vec                 lu, fieldVals;
264:   PetscScalar        *vals;
265:   const PetscScalar  *ivals, *vcoords;
266:   PetscReal          *pcoords;
267:   PetscBool           pointsAllProcs = PETSC_TRUE;
268:   PetscInt            dim, spaceDim, Nc, c, Np, p;
269:   PetscMPIInt         rank, size;
270:   PetscViewer         selfviewer;

272:   PetscFunctionBeginUser;
273:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
274:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &ctx));
275:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &ctx, &dm));
276:   PetscCall(DMGetDimension(dm, &dim));
277:   PetscCall(DMGetCoordinateDim(dm, &spaceDim));
278:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
279:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
280:   /* Create points */
281:   PetscCall(CreatePoints(dm, &Np, &pcoords, &pointsAllProcs, &ctx));
282:   /* Create interpolator */
283:   PetscCall(DMInterpolationCreate(PETSC_COMM_WORLD, &interpolator));
284:   PetscCall(DMInterpolationSetDim(interpolator, spaceDim));
285:   PetscCall(DMInterpolationAddPoints(interpolator, Np, pcoords));
286:   PetscCall(DMInterpolationSetUp(interpolator, dm, pointsAllProcs, PETSC_FALSE));
287:   /* Check locations */
288:   for (c = 0; c < interpolator->n; ++c) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %" PetscInt_FMT " is in Cell %" PetscInt_FMT "\n", rank, c, interpolator->cells[c]));
289:   PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL));
290:   PetscCall(VecView(interpolator->coords, PETSC_VIEWER_STDOUT_WORLD));
291:   Nc = dim;
292:   PetscCall(CreateDiscretization(dm, Nc, &ctx));
293:   /* Create function */
294:   PetscCall(PetscCalloc2(Nc, &funcs, Nc, &vals));
295:   switch (ctx.funcType) {
296:   case CONSTANT:
297:     for (c = 0; c < Nc; ++c) funcs[c] = constant;
298:     break;
299:   case LINEAR:
300:     for (c = 0; c < Nc; ++c) funcs[c] = linear;
301:     break;
302:   default:
303:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid function type: %d", (int)ctx.funcType);
304:   }
305:   PetscCall(DMGetLocalVector(dm, &lu));
306:   PetscCall(DMProjectFunctionLocal(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, lu));
307:   PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
308:   PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &selfviewer));
309:   PetscCall(PetscViewerASCIIPrintf(selfviewer, "[%d]solution\n", rank));
310:   PetscCall(VecView(lu, selfviewer));
311:   PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &selfviewer));
312:   PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));
313:   /* Check interpolant */
314:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, interpolator->n * Nc, &fieldVals));
315:   PetscCall(DMInterpolationSetDof(interpolator, Nc));
316:   PetscCall(DMInterpolationEvaluate(interpolator, dm, lu, fieldVals));
317:   for (p = 0; p < size; ++p) {
318:     if (p == rank) {
319:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d]Field values\n", rank));
320:       PetscCall(VecView(fieldVals, PETSC_VIEWER_STDOUT_SELF));
321:     }
322:     PetscCall(PetscBarrier((PetscObject)dm));
323:   }
324:   PetscCall(VecGetArrayRead(interpolator->coords, &vcoords));
325:   PetscCall(VecGetArrayRead(fieldVals, &ivals));
326:   for (p = 0; p < interpolator->n; ++p) {
327:     for (c = 0; c < Nc; ++c) {
328: #if defined(PETSC_USE_COMPLEX)
329:       PetscReal vcoordsReal[3];

331:       for (PetscInt i = 0; i < spaceDim; i++) vcoordsReal[i] = PetscRealPart(vcoords[p * spaceDim + i]);
332: #else
333:       const PetscReal *vcoordsReal = &vcoords[p * spaceDim];
334: #endif
335:       PetscCall((*funcs[c])(dim, 0.0, vcoordsReal, Nc, vals, NULL));
336:       PetscCheck(PetscAbsScalar(ivals[p * Nc + c] - vals[c]) <= PETSC_SQRT_MACHINE_EPSILON, 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);
337:     }
338:   }
339:   PetscCall(VecRestoreArrayRead(interpolator->coords, &vcoords));
340:   PetscCall(VecRestoreArrayRead(fieldVals, &ivals));
341:   /* Cleanup */
342:   PetscCall(PetscFree(pcoords));
343:   PetscCall(PetscFree2(funcs, vals));
344:   PetscCall(VecDestroy(&fieldVals));
345:   PetscCall(DMRestoreLocalVector(dm, &lu));
346:   PetscCall(DMInterpolationDestroy(&interpolator));
347:   PetscCall(DMDestroy(&dm));
348:   PetscCall(PetscFinalize());
349:   return 0;
350: }

352: /*TEST

354:   testset:
355:     requires: ctetgen
356:     args: -dm_plex_dim 3 -petscspace_degree 1

358:     test:
359:       suffix: 0
360:     test:
361:       suffix: 1
362:       args: -dm_refine 2
363:     test:
364:       suffix: 2
365:       nsize: 2
366:       args: -petscpartitioner_type simple
367:     test:
368:       suffix: 3
369:       nsize: 2
370:       args: -dm_refine 2 -petscpartitioner_type simple
371:     test:
372:       suffix: 4
373:       nsize: 5
374:       args: -petscpartitioner_type simple
375:     test:
376:       suffix: 5
377:       nsize: 5
378:       args: -dm_refine 2 -petscpartitioner_type simple
379:     test:
380:       suffix: 6
381:       args: -point_type grid
382:     test:
383:       suffix: 7
384:       args: -dm_refine 2 -point_type grid
385:     test:
386:       suffix: 8
387:       nsize: 2
388:       args: -petscpartitioner_type simple -point_type grid
389:     test:
390:       suffix: 9
391:       args: -point_type grid_replicated
392:     test:
393:       suffix: 10
394:       nsize: 2
395:       args: -petscpartitioner_type simple -point_type grid_replicated
396:     test:
397:       suffix: 11
398:       nsize: 2
399:       args: -dm_refine 2 -petscpartitioner_type simple -point_type grid_replicated

401:   test:
402:     suffix: fv_0
403:     requires: triangle
404:     args: -use_fv -func_type constant

406: TEST*/