Actual source code: ex42.c

  1: static const char help[] = "Simple libCEED test to calculate surface area using 1^T M 1";

  3: /*
  4:   This is a recreation of libCeed Example 2: https://libceed.readthedocs.io/en/latest/examples/ceed/
  5: */

  7: #include <petscdmceed.h>
  8: #include <petscdmplexceed.h>
  9: #include <petscfeceed.h>
 10: #include <petscdmplex.h>
 11: #include <petscds.h>

 13: typedef struct {
 14:   PetscReal         areaExact;
 15:   CeedQFunctionUser setupgeo, apply;
 16:   const char       *setupgeofname, *applyfname;
 17: } AppCtx;

 19: typedef struct {
 20:   CeedQFunction qf_apply;
 21:   CeedOperator  op_apply;
 22:   CeedVector    qdata, uceed, vceed;
 23: } CeedData;

 25: static PetscErrorCode CeedDataDestroy(CeedData *data)
 26: {
 28:   CeedVectorDestroy(&data->qdata);
 29:   CeedVectorDestroy(&data->uceed);
 30:   CeedVectorDestroy(&data->vceed);
 31:   CeedQFunctionDestroy(&data->qf_apply);
 32:   CeedOperatorDestroy(&data->op_apply);
 33:   return 0;
 34: }

 36: CEED_QFUNCTION(Mass)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out)
 37: {
 38:   const CeedScalar *u = in[0], *qdata = in[1];
 39:   CeedScalar       *v = out[0];

 41:   CeedPragmaSIMD for (CeedInt i = 0; i < Q; ++i) v[i] = qdata[i] * u[i];

 43:   return 0;
 44: }

 46: /*
 47: // Reference (parent) 2D coordinates: X \in [-1, 1]^2
 48: //
 49: // Global physical coordinates given by the mesh (3D): xx \in [-l, l]^3
 50: //
 51: // Local physical coordinates on the manifold (2D): x \in [-l, l]^2
 52: //
 53: // Change of coordinates matrix computed by the library:
 54: //   (physical 3D coords relative to reference 2D coords)
 55: //   dxx_j/dX_i (indicial notation) [3 * 2]
 56: //
 57: // Change of coordinates x (physical 2D) relative to xx (phyisical 3D):
 58: //   dx_i/dxx_j (indicial notation) [2 * 3]
 59: //
 60: // Change of coordinates x (physical 2D) relative to X (reference 2D):
 61: //   (by chain rule)
 62: //   dx_i/dX_j = dx_i/dxx_k * dxx_k/dX_j
 63: //
 64: // The quadrature data is stored in the array qdata.
 65: //
 66: // We require the determinant of the Jacobian to properly compute integrals of the form: int(u v)
 67: //
 68: // Qdata: w * det(dx_i/dX_j)
 69: */
 70: CEED_QFUNCTION(SetupMassGeoCube)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out)
 71: {
 72:   const CeedScalar *J = in[1], *w = in[2];
 73:   CeedScalar       *qdata = out[0];

 75:   CeedPragmaSIMD for (CeedInt i = 0; i < Q; ++i)
 76:   {
 77:     // Read dxxdX Jacobian entries, stored as [[0 3], [1 4], [2 5]]
 78:     const CeedScalar dxxdX[3][2] = {
 79:       {J[i + Q * 0], J[i + Q * 3]},
 80:       {J[i + Q * 1], J[i + Q * 4]},
 81:       {J[i + Q * 2], J[i + Q * 5]}
 82:     };
 83:     // Modulus of dxxdX column vectors
 84:     const CeedScalar modg1 = PetscSqrtReal(dxxdX[0][0] * dxxdX[0][0] + dxxdX[1][0] * dxxdX[1][0] + dxxdX[2][0] * dxxdX[2][0]);
 85:     const CeedScalar modg2 = PetscSqrtReal(dxxdX[0][1] * dxxdX[0][1] + dxxdX[1][1] * dxxdX[1][1] + dxxdX[2][1] * dxxdX[2][1]);
 86:     // Use normalized column vectors of dxxdX as rows of dxdxx
 87:     const CeedScalar dxdxx[2][3] = {
 88:       {dxxdX[0][0] / modg1, dxxdX[1][0] / modg1, dxxdX[2][0] / modg1},
 89:       {dxxdX[0][1] / modg2, dxxdX[1][1] / modg2, dxxdX[2][1] / modg2}
 90:     };

 92:     CeedScalar dxdX[2][2];
 93:     for (int j = 0; j < 2; ++j)
 94:       for (int k = 0; k < 2; ++k) {
 95:         dxdX[j][k] = 0;
 96:         for (int l = 0; l < 3; ++l) dxdX[j][k] += dxdxx[j][l] * dxxdX[l][k];
 97:       }
 98:     qdata[i + Q * 0] = (dxdX[0][0] * dxdX[1][1] - dxdX[1][0] * dxdX[0][1]) * w[i]; /* det J * weight */
 99:   }
100:   return 0;
101: }

103: /*
104: // Reference (parent) 2D coordinates: X \in [-1, 1]^2
105: //
106: // Global 3D physical coordinates given by the mesh: xx \in [-R, R]^3
107: //   with R radius of the sphere
108: //
109: // Local 3D physical coordinates on the 2D manifold: x \in [-l, l]^3
110: //   with l half edge of the cube inscribed in the sphere
111: //
112: // Change of coordinates matrix computed by the library:
113: //   (physical 3D coords relative to reference 2D coords)
114: //   dxx_j/dX_i (indicial notation) [3 * 2]
115: //
116: // Change of coordinates x (on the 2D manifold) relative to xx (phyisical 3D):
117: //   dx_i/dxx_j (indicial notation) [3 * 3]
118: //
119: // Change of coordinates x (on the 2D manifold) relative to X (reference 2D):
120: //   (by chain rule)
121: //   dx_i/dX_j = dx_i/dxx_k * dxx_k/dX_j [3 * 2]
122: //
123: // modJ is given by the magnitude of the cross product of the columns of dx_i/dX_j
124: //
125: // The quadrature data is stored in the array qdata.
126: //
127: // We require the determinant of the Jacobian to properly compute integrals of
128: //   the form: int(u v)
129: //
130: // Qdata: modJ * w
131: */
132: CEED_QFUNCTION(SetupMassGeoSphere)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out)
133: {
134:   const CeedScalar *X = in[0], *J = in[1], *w = in[2];
135:   CeedScalar       *qdata = out[0];

137:   CeedPragmaSIMD for (CeedInt i = 0; i < Q; ++i)
138:   {
139:     const CeedScalar xx[3][1] = {{X[i + 0 * Q]}, {X[i + 1 * Q]}, {X[i + 2 * Q]}};
140:     // Read dxxdX Jacobian entries, stored as [[0 3], [1 4], [2 5]]
141:     const CeedScalar dxxdX[3][2] = {
142:       {J[i + Q * 0], J[i + Q * 3]},
143:       {J[i + Q * 1], J[i + Q * 4]},
144:       {J[i + Q * 2], J[i + Q * 5]}
145:     };
146:     // Setup
147:     const CeedScalar modxxsq = xx[0][0] * xx[0][0] + xx[1][0] * xx[1][0] + xx[2][0] * xx[2][0];
148:     CeedScalar       xxsq[3][3];
149:     for (int j = 0; j < 3; ++j)
150:       for (int k = 0; k < 3; ++k) {
151:         xxsq[j][k] = 0.;
152:         for (int l = 0; l < 1; ++l) xxsq[j][k] += xx[j][l] * xx[k][l] / (sqrt(modxxsq) * modxxsq);
153:       }

155:     const CeedScalar dxdxx[3][3] = {
156:       {1. / sqrt(modxxsq) - xxsq[0][0], -xxsq[0][1],                     -xxsq[0][2]                    },
157:       {-xxsq[1][0],                     1. / sqrt(modxxsq) - xxsq[1][1], -xxsq[1][2]                    },
158:       {-xxsq[2][0],                     -xxsq[2][1],                     1. / sqrt(modxxsq) - xxsq[2][2]}
159:     };

161:     CeedScalar dxdX[3][2];
162:     for (int j = 0; j < 3; ++j)
163:       for (int k = 0; k < 2; ++k) {
164:         dxdX[j][k] = 0.;
165:         for (int l = 0; l < 3; ++l) dxdX[j][k] += dxdxx[j][l] * dxxdX[l][k];
166:       }
167:     // J is given by the cross product of the columns of dxdX
168:     const CeedScalar J[3][1] = {{dxdX[1][0] * dxdX[2][1] - dxdX[2][0] * dxdX[1][1]}, {dxdX[2][0] * dxdX[0][1] - dxdX[0][0] * dxdX[2][1]}, {dxdX[0][0] * dxdX[1][1] - dxdX[1][0] * dxdX[0][1]}};
169:     // Use the magnitude of J as our detJ (volume scaling factor)
170:     const CeedScalar modJ = sqrt(J[0][0] * J[0][0] + J[1][0] * J[1][0] + J[2][0] * J[2][0]);
171:     qdata[i + Q * 0]      = modJ * w[i];
172:   }
173:   return 0;
174: }

176: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *ctx)
177: {
178:   DMPlexShape shape = DM_SHAPE_UNKNOWN;

181:   PetscOptionsBegin(comm, "", "libCEED Test Options", "DMPLEX");
182:   PetscOptionsEnd();
183:   PetscOptionsGetEnum(NULL, NULL, "-dm_plex_shape", DMPlexShapes, (PetscEnum *)&shape, NULL);
184:   ctx->setupgeo      = NULL;
185:   ctx->setupgeofname = NULL;
186:   ctx->apply         = Mass;
187:   ctx->applyfname    = Mass_loc;
188:   ctx->areaExact     = 0.0;
189:   switch (shape) {
190:   case DM_SHAPE_BOX_SURFACE:
191:     ctx->setupgeo      = SetupMassGeoCube;
192:     ctx->setupgeofname = SetupMassGeoCube_loc;
193:     ctx->areaExact     = 6.0;
194:     break;
195:   case DM_SHAPE_SPHERE:
196:     ctx->setupgeo      = SetupMassGeoSphere;
197:     ctx->setupgeofname = SetupMassGeoSphere_loc;
198:     ctx->areaExact     = 4.0 * M_PI;
199:     break;
200:   default:
201:     break;
202:   }
203:   return 0;
204: }

206: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm)
207: {
208:   DMCreate(comm, dm);
209:   DMSetType(*dm, DMPLEX);
210:   DMSetFromOptions(*dm);
211:   DMViewFromOptions(*dm, NULL, "-dm_view");
212: #ifdef PETSC_HAVE_LIBCEED
213:   {
214:     Ceed        ceed;
215:     const char *usedresource;

217:     DMGetCeed(*dm, &ceed);
218:     CeedGetResource(ceed, &usedresource);
219:     PetscPrintf(PetscObjectComm((PetscObject)*dm), "libCEED Backend: %s\n", usedresource);
220:   }
221: #endif
222:   return 0;
223: }

225: static PetscErrorCode SetupDiscretization(DM dm)
226: {
227:   DM        cdm;
228:   PetscFE   fe, cfe;
229:   PetscInt  dim, cnc;
230:   PetscBool simplex;

233:   DMGetDimension(dm, &dim);
234:   DMPlexIsSimplex(dm, &simplex);
235:   PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, NULL, PETSC_DETERMINE, &fe);
236:   PetscFESetName(fe, "indicator");
237:   DMAddField(dm, NULL, (PetscObject)fe);
238:   PetscFEDestroy(&fe);
239:   DMCreateDS(dm);
240:   DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL);
241:   DMGetCoordinateDim(dm, &cnc);
242:   PetscFECreateDefault(PETSC_COMM_SELF, dim, cnc, simplex, NULL, PETSC_DETERMINE, &cfe);
243:   DMProjectCoordinates(dm, cfe);
244:   PetscFEDestroy(&cfe);
245:   DMGetCoordinateDM(dm, &cdm);
246:   DMPlexSetClosurePermutationTensor(cdm, PETSC_DETERMINE, NULL);
247:   return 0;
248: }

250: static PetscErrorCode LibCeedSetupByDegree(DM dm, AppCtx *ctx, CeedData *data)
251: {
252:   PetscDS             ds;
253:   PetscFE             fe, cfe;
254:   Ceed                ceed;
255:   CeedElemRestriction Erestrictx, Erestrictu, Erestrictq;
256:   CeedQFunction       qf_setupgeo;
257:   CeedOperator        op_setupgeo;
258:   CeedVector          xcoord;
259:   CeedBasis           basisu, basisx;
260:   CeedInt             Nqdata = 1;
261:   CeedInt             nqpts, nqptsx;
262:   DM                  cdm;
263:   Vec                 coords;
264:   const PetscScalar  *coordArray;
265:   PetscInt            dim, cdim, cStart, cEnd, Ncell;

268:   DMGetCeed(dm, &ceed);
269:   DMGetDimension(dm, &dim);
270:   DMGetCoordinateDim(dm, &cdim);
271:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
272:   Ncell = cEnd - cStart;
273:   // CEED bases
274:   DMGetDS(dm, &ds);
275:   PetscDSGetDiscretization(ds, 0, (PetscObject *)&fe);
276:   PetscFEGetCeedBasis(fe, &basisu);
277:   DMGetCoordinateDM(dm, &cdm);
278:   DMGetDS(cdm, &ds);
279:   PetscDSGetDiscretization(ds, 0, (PetscObject *)&cfe);
280:   PetscFEGetCeedBasis(cfe, &basisx);

282:   DMPlexGetCeedRestriction(cdm, NULL, 0, 0, 0, &Erestrictx);
283:   DMPlexGetCeedRestriction(dm, NULL, 0, 0, 0, &Erestrictu);
284:   CeedBasisGetNumQuadraturePoints(basisu, &nqpts);
285:   CeedBasisGetNumQuadraturePoints(basisx, &nqptsx);
287:   CeedElemRestrictionCreateStrided(ceed, Ncell, nqpts, Nqdata, Nqdata * Ncell * nqpts, CEED_STRIDES_BACKEND, &Erestrictq);

289:   DMGetCoordinatesLocal(dm, &coords);
290:   VecGetArrayRead(coords, &coordArray);
291:   CeedElemRestrictionCreateVector(Erestrictx, &xcoord, NULL);
292:   CeedVectorSetArray(xcoord, CEED_MEM_HOST, CEED_COPY_VALUES, (PetscScalar *)coordArray);
293:   VecRestoreArrayRead(coords, &coordArray);

295:   // Create the vectors that will be needed in setup and apply
296:   CeedElemRestrictionCreateVector(Erestrictu, &data->uceed, NULL);
297:   CeedElemRestrictionCreateVector(Erestrictu, &data->vceed, NULL);
298:   CeedElemRestrictionCreateVector(Erestrictq, &data->qdata, NULL);

300:   // Create the Q-function that builds the operator (i.e. computes its quadrature data) and set its context data
301:   CeedQFunctionCreateInterior(ceed, 1, ctx->setupgeo, ctx->setupgeofname, &qf_setupgeo);
302:   CeedQFunctionAddInput(qf_setupgeo, "x", cdim, CEED_EVAL_INTERP);
303:   CeedQFunctionAddInput(qf_setupgeo, "dx", cdim * dim, CEED_EVAL_GRAD);
304:   CeedQFunctionAddInput(qf_setupgeo, "weight", 1, CEED_EVAL_WEIGHT);
305:   CeedQFunctionAddOutput(qf_setupgeo, "qdata", Nqdata, CEED_EVAL_NONE);

307:   // Set up the mass operator
308:   CeedQFunctionCreateInterior(ceed, 1, ctx->apply, ctx->applyfname, &data->qf_apply);
309:   CeedQFunctionAddInput(data->qf_apply, "u", 1, CEED_EVAL_INTERP);
310:   CeedQFunctionAddInput(data->qf_apply, "qdata", Nqdata, CEED_EVAL_NONE);
311:   CeedQFunctionAddOutput(data->qf_apply, "v", 1, CEED_EVAL_INTERP);

313:   // Create the operator that builds the quadrature data for the operator
314:   CeedOperatorCreate(ceed, qf_setupgeo, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setupgeo);
315:   CeedOperatorSetField(op_setupgeo, "x", Erestrictx, basisx, CEED_VECTOR_ACTIVE);
316:   CeedOperatorSetField(op_setupgeo, "dx", Erestrictx, basisx, CEED_VECTOR_ACTIVE);
317:   CeedOperatorSetField(op_setupgeo, "weight", CEED_ELEMRESTRICTION_NONE, basisx, CEED_VECTOR_NONE);
318:   CeedOperatorSetField(op_setupgeo, "qdata", Erestrictq, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);

320:   // Create the mass operator
321:   CeedOperatorCreate(ceed, data->qf_apply, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &data->op_apply);
322:   CeedOperatorSetField(data->op_apply, "u", Erestrictu, basisu, CEED_VECTOR_ACTIVE);
323:   CeedOperatorSetField(data->op_apply, "qdata", Erestrictq, CEED_BASIS_COLLOCATED, data->qdata);
324:   CeedOperatorSetField(data->op_apply, "v", Erestrictu, basisu, CEED_VECTOR_ACTIVE);

326:   // Setup qdata
327:   CeedOperatorApply(op_setupgeo, xcoord, data->qdata, CEED_REQUEST_IMMEDIATE);

329:   CeedElemRestrictionDestroy(&Erestrictq);
330:   CeedQFunctionDestroy(&qf_setupgeo);
331:   CeedOperatorDestroy(&op_setupgeo);
332:   CeedVectorDestroy(&xcoord);
333:   return 0;
334: }

336: int main(int argc, char **argv)
337: {
338:   MPI_Comm     comm;
339:   DM           dm;
340:   AppCtx       ctx;
341:   Vec          U, Uloc, V, Vloc;
342:   PetscScalar *v;
343:   PetscScalar  area;
344:   CeedData     ceeddata;

347:   PetscInitialize(&argc, &argv, NULL, help);
348:   comm = PETSC_COMM_WORLD;
349:   ProcessOptions(comm, &ctx);
350:   CreateMesh(comm, &ctx, &dm);
351:   SetupDiscretization(dm);

353:   LibCeedSetupByDegree(dm, &ctx, &ceeddata);

355:   DMCreateGlobalVector(dm, &U);
356:   DMCreateLocalVector(dm, &Uloc);
357:   VecDuplicate(U, &V);
358:   VecDuplicate(Uloc, &Vloc);

360:   /**/
361:   VecSet(Uloc, 1.);
362:   VecZeroEntries(V);
363:   VecZeroEntries(Vloc);
364:   VecGetArray(Vloc, &v);
365:   CeedVectorSetArray(ceeddata.vceed, CEED_MEM_HOST, CEED_USE_POINTER, v);
366:   CeedVectorSetValue(ceeddata.uceed, 1.0);
367:   CeedOperatorApply(ceeddata.op_apply, ceeddata.uceed, ceeddata.vceed, CEED_REQUEST_IMMEDIATE);
368:   CeedVectorTakeArray(ceeddata.vceed, CEED_MEM_HOST, NULL);
369:   VecRestoreArray(Vloc, &v);
370:   DMLocalToGlobalBegin(dm, Vloc, ADD_VALUES, V);
371:   DMLocalToGlobalEnd(dm, Vloc, ADD_VALUES, V);

373:   VecSum(V, &area);
374:   if (ctx.areaExact > 0.) {
375:     PetscReal error = PetscAbsReal(area - ctx.areaExact);
376:     PetscReal tol   = PETSC_SMALL;

378:     PetscPrintf(comm, "Exact mesh surface area    : % .*f\n", PetscAbsReal(ctx.areaExact - round(ctx.areaExact)) > 1E-15 ? 14 : 1, (double)ctx.areaExact);
379:     PetscPrintf(comm, "Computed mesh surface area : % .*f\n", PetscAbsScalar(area - round(area)) > 1E-15 ? 14 : 1, (double)PetscRealPart(area));
380:     if (error > tol) {
381:       PetscPrintf(comm, "Area error                 : % .14g\n", (double)error);
382:     } else {
383:       PetscPrintf(comm, "Area verifies!\n");
384:     }
385:   }

387:   CeedDataDestroy(&ceeddata);
388:   VecDestroy(&U);
389:   VecDestroy(&Uloc);
390:   VecDestroy(&V);
391:   VecDestroy(&Vloc);
392:   DMDestroy(&dm);
393:   return PetscFinalize();
394: }

396: /*TEST

398:   build:
399:     requires: libceed

401:   testset:
402:     args: -dm_plex_simplex 0 -petscspace_degree 3 -dm_view -dm_petscds_view \
403:           -petscfe_default_quadrature_order 4 -coord_dm_default_quadrature_order 4

405:     test:
406:       suffix: cube_3
407:       args: -dm_plex_shape box_surface -dm_refine 2

409:     test:
410:       suffix: cube_3_p4
411:       nsize: 4
412:       args: -petscpartitioner_type simple -dm_refine_pre 1 -dm_plex_shape box_surface -dm_refine 1

414:     test:
415:       suffix: sphere_3
416:       args: -dm_plex_shape sphere -dm_refine 3

418:     test:
419:       suffix: sphere_3_p4
420:       nsize: 4
421:       args: -petscpartitioner_type simple -dm_refine_pre 1 -dm_plex_shape sphere -dm_refine 2

423: TEST*/