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: {
 27:   PetscFunctionBeginUser;
 28:   PetscCall(CeedVectorDestroy(&data->qdata));
 29:   PetscCall(CeedVectorDestroy(&data->uceed));
 30:   PetscCall(CeedVectorDestroy(&data->vceed));
 31:   PetscCall(CeedQFunctionDestroy(&data->qf_apply));
 32:   PetscCall(CeedOperatorDestroy(&data->op_apply));
 33:   PetscFunctionReturn(PETSC_SUCCESS);
 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;

180:   PetscFunctionBeginUser;
181:   PetscOptionsBegin(comm, "", "libCEED Test Options", "DMPLEX");
182:   PetscOptionsEnd();
183:   PetscCall(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:   PetscFunctionReturn(PETSC_SUCCESS);
204: }

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

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

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

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

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

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

283:   PetscCall(DMPlexGetCeedRestriction(cdm, NULL, 0, 0, 0, &Erestrictx));
284:   PetscCall(DMPlexGetCeedRestriction(dm, NULL, 0, 0, 0, &Erestrictu));
285:   PetscCall(CeedBasisGetNumQuadraturePoints(basisu, &nqpts));
286:   PetscCall(CeedBasisGetNumQuadraturePoints(basisx, &nqptsx));
287:   PetscCheck(nqptsx == nqpts, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of qpoints for u %" PetscInt_FMT " != %" PetscInt_FMT " Number of qpoints for x", nqpts, nqptsx);
288:   PetscCall(CeedElemRestrictionCreateStrided(ceed, Ncell, nqpts, Nqdata, Nqdata * Ncell * nqpts, CEED_STRIDES_BACKEND, &Erestrictq));

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

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

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

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

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

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

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

330:   PetscCall(CeedElemRestrictionDestroy(&Erestrictq));
331:   PetscCall(CeedQFunctionDestroy(&qf_setupgeo));
332:   PetscCall(CeedOperatorDestroy(&op_setupgeo));
333:   PetscCall(CeedVectorDestroy(&xcoord));
334:   PetscFunctionReturn(PETSC_SUCCESS);
335: }

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

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

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

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

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

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

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

387:   PetscCall(CeedDataDestroy(&ceeddata));
388:   PetscCall(VecDestroy(&U));
389:   PetscCall(VecDestroy(&Uloc));
390:   PetscCall(VecDestroy(&V));
391:   PetscCall(VecDestroy(&Vloc));
392:   PetscCall(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 -cdm_default_quadrature_order 4
404:     filter: sed "s /cpu/self/xsmm /cpu/self/avx "

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

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

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

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

424: TEST*/