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 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 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]; 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 */ 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]; 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 */ 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]; 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 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 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 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(DMProjectCoordinates(dm, cfe)); 245 PetscCall(PetscFEDestroy(&cfe)); 246 PetscCall(DMGetCoordinateDM(dm, &cdm)); 247 PetscCall(DMPlexSetClosurePermutationTensor(cdm, PETSC_DETERMINE, NULL)); 248 PetscFunctionReturn(PETSC_SUCCESS); 249 } 250 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 %" 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)); 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_COLLOCATED, 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_COLLOCATED, 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 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 "s /cpu/self/xsmm /cpu/self/avx " 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