1 static const char help[] = "Simple libCEED test to calculate surface area using 1^T M 1";
2
3 /*
4 This is a recreation of libCeed Example 2: https://libceed.readthedocs.io/en/latest/examples/ceed/
5 */
6
7 #include <petscdmceed.h>
8 #include <petscdmplexceed.h>
9 #include <petscfeceed.h>
10 #include <petscdmplex.h>
11 #include <petscds.h>
12
13 typedef struct {
14 PetscReal areaExact;
15 CeedQFunctionUser setupgeo, apply;
16 const char *setupgeofname, *applyfname;
17 } AppCtx;
18
19 typedef struct {
20 CeedQFunction qf_apply;
21 CeedOperator op_apply;
22 CeedVector qdata, uceed, vceed;
23 } CeedData;
24
CeedDataDestroy(CeedData * data)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 }
35
Mass(PetscCtx ctx,const CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)36 CEED_QFUNCTION(Mass)(PetscCtx 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];
40
41 CeedPragmaSIMD for (CeedInt i = 0; i < Q; ++i) v[i] = qdata[i] * u[i];
42
43 return 0;
44 }
45
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 */
SetupMassGeoCube(PetscCtx ctx,const CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)70 CEED_QFUNCTION(SetupMassGeoCube)(PetscCtx 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];
74
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 };
91
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 }
102
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 */
SetupMassGeoSphere(PetscCtx ctx,const CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)132 CEED_QFUNCTION(SetupMassGeoSphere)(PetscCtx 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];
136
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 }
154
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 };
160
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 }
175
ProcessOptions(MPI_Comm comm,AppCtx * ctx)176 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *ctx)
177 {
178 DMPlexShape shape = DM_SHAPE_UNKNOWN;
179
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 }
205
CreateMesh(MPI_Comm comm,AppCtx * ctx,DM * dm)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;
217
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 }
225
SetupDiscretization(DM dm)226 static PetscErrorCode SetupDiscretization(DM dm)
227 {
228 DM cdm;
229 PetscFE fe, cfe;
230 PetscInt dim, cnc;
231 PetscBool simplex;
232
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(DMSetCoordinateDisc(dm, cfe, PETSC_FALSE, PETSC_TRUE));
245 PetscCall(PetscFEDestroy(&cfe));
246 PetscCall(DMGetCoordinateDM(dm, &cdm));
247 PetscCall(DMPlexSetClosurePermutationTensor(cdm, PETSC_DETERMINE, NULL));
248 PetscFunctionReturn(PETSC_SUCCESS);
249 }
250
LibCeedSetupByDegree(DM dm,AppCtx * ctx,CeedData * data)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;
267
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));
282
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 %" CeedInt_FMT " != %" CeedInt_FMT " Number of qpoints for x", nqpts, nqptsx);
288 PetscCall(CeedElemRestrictionCreateStrided(ceed, Ncell, nqpts, Nqdata, Nqdata * Ncell * nqpts, CEED_STRIDES_BACKEND, &Erestrictq));
289
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));
295
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));
300
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));
307
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));
313
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_NONE, CEED_VECTOR_ACTIVE));
320
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_NONE, data->qdata));
325 PetscCall(CeedOperatorSetField(data->op_apply, "v", Erestrictu, basisu, CEED_VECTOR_ACTIVE));
326
327 // Setup qdata
328 PetscCall(CeedOperatorApply(op_setupgeo, xcoord, data->qdata, CEED_REQUEST_IMMEDIATE));
329
330 PetscCall(CeedElemRestrictionDestroy(&Erestrictq));
331 PetscCall(CeedQFunctionDestroy(&qf_setupgeo));
332 PetscCall(CeedOperatorDestroy(&op_setupgeo));
333 PetscCall(CeedVectorDestroy(&xcoord));
334 PetscFunctionReturn(PETSC_SUCCESS);
335 }
336
main(int argc,char ** argv)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;
346
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));
352
353 PetscCall(LibCeedSetupByDegree(dm, &ctx, &ceeddata));
354
355 PetscCall(DMCreateGlobalVector(dm, &U));
356 PetscCall(DMCreateLocalVector(dm, &Uloc));
357 PetscCall(VecDuplicate(U, &V));
358 PetscCall(VecDuplicate(Uloc, &Vloc));
359
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));
372
373 PetscCall(VecSum(V, &area));
374 if (ctx.areaExact > 0.) {
375 PetscReal error = PetscAbsReal(area - ctx.areaExact);
376 PetscReal tol = PETSC_SMALL;
377
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 }
386
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 }
395
396 /*TEST
397
398 build:
399 requires: libceed
400
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 -e "s /cpu/self/xsmm /cpu/self/opt " -e "s /cpu/self/avx /cpu/self/opt "
405
406 test:
407 suffix: cube_3
408 args: -dm_plex_shape box_surface -dm_refine 2
409
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
414
415 test:
416 suffix: sphere_3
417 args: -dm_plex_shape sphere -dm_refine 3
418
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
423
424 TEST*/
425