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*/