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 struct {
 15:   PointType pointType; /* Point generation mechanism */
 16: } AppCtx;

 18: static PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 19: {
 20:   PetscInt d, c;

 22:   PetscFunctionBeginUser;
 23:   PetscCheck(Nc == 3, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Something is wrong: %" PetscInt_FMT, Nc);
 24:   for (c = 0; c < Nc; ++c) {
 25:     u[c] = 0.0;
 26:     for (d = 0; d < dim; ++d) u[c] += x[d];
 27:   }
 28:   PetscFunctionReturn(PETSC_SUCCESS);
 29: }

 31: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 32: {
 33:   const char *pointTypes[3] = {"centroid", "grid", "grid_replicated"};
 34:   PetscInt    pt;

 36:   PetscFunctionBegin;
 37:   options->pointType = CENTROID;
 38:   PetscOptionsBegin(comm, "", "Interpolation Options", "DMPLEX");
 39:   pt = options->pointType;
 40:   PetscCall(PetscOptionsEList("-point_type", "The point type", "ex2.c", pointTypes, 3, pointTypes[options->pointType], &pt, NULL));
 41:   options->pointType = (PointType)pt;
 42:   PetscOptionsEnd();
 43:   PetscFunctionReturn(PETSC_SUCCESS);
 44: }

 46: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm)
 47: {
 48:   PetscFunctionBegin;
 49:   PetscCall(DMCreate(comm, dm));
 50:   PetscCall(DMSetType(*dm, DMPLEX));
 51:   PetscCall(DMSetFromOptions(*dm));
 52:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
 53:   PetscFunctionReturn(PETSC_SUCCESS);
 54: }

 56: static PetscErrorCode CreatePoints_Centroid(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx)
 57: {
 58:   PetscSection coordSection;
 59:   Vec          coordsLocal;
 60:   PetscInt     spaceDim, p;
 61:   PetscMPIInt  rank;

 63:   PetscFunctionBegin;
 64:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 65:   PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
 66:   PetscCall(DMGetCoordinateSection(dm, &coordSection));
 67:   PetscCall(DMGetCoordinateDim(dm, &spaceDim));
 68:   PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, Np));
 69:   PetscCall(PetscCalloc1(*Np * spaceDim, pcoords));
 70:   for (p = 0; p < *Np; ++p) {
 71:     PetscScalar *coords = NULL;
 72:     PetscInt     size, num, n, d;

 74:     PetscCall(DMPlexVecGetClosure(dm, coordSection, coordsLocal, p, &size, &coords));
 75:     num = size / spaceDim;
 76:     for (n = 0; n < num; ++n) {
 77:       for (d = 0; d < spaceDim; ++d) (*pcoords)[p * spaceDim + d] += PetscRealPart(coords[n * spaceDim + d]) / num;
 78:     }
 79:     PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %" PetscInt_FMT " (", rank, p));
 80:     for (d = 0; d < spaceDim; ++d) {
 81:       PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%g", (double)(*pcoords)[p * spaceDim + d]));
 82:       if (d < spaceDim - 1) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ", "));
 83:     }
 84:     PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ")\n"));
 85:     PetscCall(DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, p, &num, &coords));
 86:   }
 87:   PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL));
 88:   *pointsAllProcs = PETSC_FALSE;
 89:   PetscFunctionReturn(PETSC_SUCCESS);
 90: }

 92: static PetscErrorCode CreatePoints_Grid(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx)
 93: {
 94:   DM            da;
 95:   DMDALocalInfo info;
 96:   PetscInt      N = 3, n = 0, dim, spaceDim, i, j, k, *ind, d;
 97:   PetscReal    *h;
 98:   PetscMPIInt   rank;

100:   PetscFunctionBegin;
101:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
102:   PetscCall(DMGetDimension(dm, &dim));
103:   PetscCall(DMGetCoordinateDim(dm, &spaceDim));
104:   PetscCall(PetscCalloc1(spaceDim, &ind));
105:   PetscCall(PetscCalloc1(spaceDim, &h));
106:   h[0] = 1.0 / (N - 1);
107:   h[1] = 1.0 / (N - 1);
108:   h[2] = 1.0 / (N - 1);
109:   PetscCall(DMDACreate(PetscObjectComm((PetscObject)dm), &da));
110:   PetscCall(DMSetDimension(da, dim));
111:   PetscCall(DMDASetSizes(da, N, N, N));
112:   PetscCall(DMDASetDof(da, 1));
113:   PetscCall(DMDASetStencilWidth(da, 1));
114:   PetscCall(DMSetUp(da));
115:   PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
116:   PetscCall(DMDAGetLocalInfo(da, &info));
117:   *Np = info.xm * info.ym * info.zm;
118:   PetscCall(PetscCalloc1(*Np * spaceDim, pcoords));
119:   for (k = info.zs; k < info.zs + info.zm; ++k) {
120:     ind[2] = k;
121:     for (j = info.ys; j < info.ys + info.ym; ++j) {
122:       ind[1] = j;
123:       for (i = info.xs; i < info.xs + info.xm; ++i, ++n) {
124:         ind[0] = i;

126:         for (d = 0; d < spaceDim; ++d) (*pcoords)[n * spaceDim + d] = ind[d] * h[d];
127:         PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %" PetscInt_FMT " (", rank, n));
128:         for (d = 0; d < spaceDim; ++d) {
129:           PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%g", (double)(*pcoords)[n * spaceDim + d]));
130:           if (d < spaceDim - 1) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ", "));
131:         }
132:         PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ")\n"));
133:       }
134:     }
135:   }
136:   PetscCall(DMDestroy(&da));
137:   PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL));
138:   PetscCall(PetscFree(ind));
139:   PetscCall(PetscFree(h));
140:   *pointsAllProcs = PETSC_FALSE;
141:   PetscFunctionReturn(PETSC_SUCCESS);
142: }

144: static PetscErrorCode CreatePoints_GridReplicated(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx)
145: {
146:   PetscInt    N = 3, n = 0, dim, spaceDim, i, j, k, *ind, d;
147:   PetscReal  *h;
148:   PetscMPIInt rank;

150:   PetscFunctionBeginUser;
151:   PetscCall(DMGetDimension(dm, &dim));
152:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
153:   PetscCall(DMGetCoordinateDim(dm, &spaceDim));
154:   PetscCall(PetscCalloc1(spaceDim, &ind));
155:   PetscCall(PetscCalloc1(spaceDim, &h));
156:   h[0] = 1.0 / (N - 1);
157:   h[1] = 1.0 / (N - 1);
158:   h[2] = 1.0 / (N - 1);
159:   *Np  = N * (dim > 1 ? N : 1) * (dim > 2 ? N : 1);
160:   PetscCall(PetscCalloc1(*Np * spaceDim, pcoords));
161:   for (k = 0; k < N; ++k) {
162:     ind[2] = k;
163:     for (j = 0; j < N; ++j) {
164:       ind[1] = j;
165:       for (i = 0; i < N; ++i, ++n) {
166:         ind[0] = i;

168:         for (d = 0; d < spaceDim; ++d) (*pcoords)[n * spaceDim + d] = ind[d] * h[d];
169:         PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %" PetscInt_FMT " (", rank, n));
170:         for (d = 0; d < spaceDim; ++d) {
171:           PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%g", (double)(*pcoords)[n * spaceDim + d]));
172:           if (d < spaceDim - 1) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ", "));
173:         }
174:         PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, ")\n"));
175:       }
176:     }
177:   }
178:   PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL));
179:   *pointsAllProcs = PETSC_TRUE;
180:   PetscCall(PetscFree(ind));
181:   PetscCall(PetscFree(h));
182:   PetscFunctionReturn(PETSC_SUCCESS);
183: }

185: static PetscErrorCode CreatePoints(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx)
186: {
187:   PetscFunctionBegin;
188:   *pointsAllProcs = PETSC_FALSE;
189:   switch (ctx->pointType) {
190:   case CENTROID:
191:     PetscCall(CreatePoints_Centroid(dm, Np, pcoords, pointsAllProcs, ctx));
192:     break;
193:   case GRID:
194:     PetscCall(CreatePoints_Grid(dm, Np, pcoords, pointsAllProcs, ctx));
195:     break;
196:   case GRID_REPLICATED:
197:     PetscCall(CreatePoints_GridReplicated(dm, Np, pcoords, pointsAllProcs, ctx));
198:     break;
199:   default:
200:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid point generation type %d", (int)ctx->pointType);
201:   }
202:   PetscFunctionReturn(PETSC_SUCCESS);
203: }

205: int main(int argc, char **argv)
206: {
207:   AppCtx ctx;
208:   PetscErrorCode (**funcs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
209:   DM                  dm;
210:   PetscFE             fe;
211:   DMInterpolationInfo interpolator;
212:   Vec                 lu, fieldVals;
213:   PetscScalar        *vals;
214:   const PetscScalar  *ivals, *vcoords;
215:   PetscReal          *pcoords;
216:   PetscBool           simplex, pointsAllProcs = PETSC_TRUE;
217:   PetscInt            dim, spaceDim, Nc, c, Np, p;
218:   PetscMPIInt         rank, size;
219:   PetscViewer         selfviewer;

221:   PetscFunctionBeginUser;
222:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
223:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &ctx));
224:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &ctx, &dm));
225:   PetscCall(DMGetDimension(dm, &dim));
226:   PetscCall(DMGetCoordinateDim(dm, &spaceDim));
227:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
228:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
229:   /* Create points */
230:   PetscCall(CreatePoints(dm, &Np, &pcoords, &pointsAllProcs, &ctx));
231:   /* Create interpolator */
232:   PetscCall(DMInterpolationCreate(PETSC_COMM_WORLD, &interpolator));
233:   PetscCall(DMInterpolationSetDim(interpolator, spaceDim));
234:   PetscCall(DMInterpolationAddPoints(interpolator, Np, pcoords));
235:   PetscCall(DMInterpolationSetUp(interpolator, dm, pointsAllProcs, PETSC_FALSE));
236:   /* Check locations */
237:   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]));
238:   PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL));
239:   PetscCall(VecView(interpolator->coords, PETSC_VIEWER_STDOUT_WORLD));
240:   /* Setup Discretization */
241:   Nc = dim;
242:   PetscCall(DMPlexIsSimplex(dm, &simplex));
243:   PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, Nc, simplex, NULL, -1, &fe));
244:   PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
245:   PetscCall(DMCreateDS(dm));
246:   PetscCall(PetscFEDestroy(&fe));
247:   /* Create function */
248:   PetscCall(PetscCalloc2(Nc, &funcs, Nc, &vals));
249:   for (c = 0; c < Nc; ++c) funcs[c] = linear;
250:   PetscCall(DMGetLocalVector(dm, &lu));
251:   PetscCall(DMProjectFunctionLocal(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, lu));
252:   PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
253:   PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &selfviewer));
254:   PetscCall(PetscViewerASCIIPrintf(selfviewer, "[%d]solution\n", rank));
255:   PetscCall(VecView(lu, selfviewer));
256:   PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &selfviewer));
257:   PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
258:   PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));
259:   /* Check interpolant */
260:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, interpolator->n * Nc, &fieldVals));
261:   PetscCall(DMInterpolationSetDof(interpolator, Nc));
262:   PetscCall(DMInterpolationEvaluate(interpolator, dm, lu, fieldVals));
263:   for (p = 0; p < size; ++p) {
264:     if (p == rank) {
265:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d]Field values\n", rank));
266:       PetscCall(VecView(fieldVals, PETSC_VIEWER_STDOUT_SELF));
267:     }
268:     PetscCall(PetscBarrier((PetscObject)dm));
269:   }
270:   PetscCall(VecGetArrayRead(interpolator->coords, &vcoords));
271:   PetscCall(VecGetArrayRead(fieldVals, &ivals));
272:   for (p = 0; p < interpolator->n; ++p) {
273:     for (c = 0; c < Nc; ++c) {
274: #if defined(PETSC_USE_COMPLEX)
275:       PetscReal vcoordsReal[3];

277:       for (PetscInt i = 0; i < spaceDim; i++) vcoordsReal[i] = PetscRealPart(vcoords[p * spaceDim + i]);
278: #else
279:       const PetscReal *vcoordsReal = &vcoords[p * spaceDim];
280: #endif
281:       PetscCall((*funcs[c])(dim, 0.0, vcoordsReal, Nc, vals, NULL));
282:       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);
283:     }
284:   }
285:   PetscCall(VecRestoreArrayRead(interpolator->coords, &vcoords));
286:   PetscCall(VecRestoreArrayRead(fieldVals, &ivals));
287:   /* Cleanup */
288:   PetscCall(PetscFree(pcoords));
289:   PetscCall(PetscFree2(funcs, vals));
290:   PetscCall(VecDestroy(&fieldVals));
291:   PetscCall(DMRestoreLocalVector(dm, &lu));
292:   PetscCall(DMInterpolationDestroy(&interpolator));
293:   PetscCall(DMDestroy(&dm));
294:   PetscCall(PetscFinalize());
295:   return 0;
296: }

298: /*TEST

300:   testset:
301:     requires: ctetgen
302:     args: -dm_plex_dim 3 -petscspace_degree 1

304:     test:
305:       suffix: 0
306:     test:
307:       suffix: 1
308:       args: -dm_refine 2
309:     test:
310:       suffix: 2
311:       nsize: 2
312:       args: -petscpartitioner_type simple
313:     test:
314:       suffix: 3
315:       nsize: 2
316:       args: -dm_refine 2 -petscpartitioner_type simple
317:     test:
318:       suffix: 4
319:       nsize: 5
320:       args: -petscpartitioner_type simple
321:     test:
322:       suffix: 5
323:       nsize: 5
324:       args: -dm_refine 2 -petscpartitioner_type simple
325:     test:
326:       suffix: 6
327:       args: -point_type grid
328:     test:
329:       suffix: 7
330:       args: -dm_refine 2 -point_type grid
331:     test:
332:       suffix: 8
333:       nsize: 2
334:       args: -petscpartitioner_type simple -point_type grid
335:     test:
336:       suffix: 9
337:       args: -point_type grid_replicated
338:     test:
339:       suffix: 10
340:       nsize: 2
341:       args: -petscpartitioner_type simple -point_type grid_replicated
342:     test:
343:       suffix: 11
344:       nsize: 2
345:       args: -dm_refine 2 -petscpartitioner_type simple -point_type grid_replicated

347: TEST*/