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