1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 // libCEED + PETSc Example: Surface Area 18 // 19 // This example demonstrates a simple usage of libCEED with PETSc to calculate 20 // the surface area of a simple closed surface, such as the one of a cube or a 21 // tensor-product discrete sphere via the mass operator. 22 // 23 // The code uses higher level communication protocols in DMPlex. 24 // 25 // Build with: 26 // 27 // make area [PETSC_DIR=</path/to/petsc>] [CEED_DIR=</path/to/libceed>] 28 // 29 // Sample runs: 30 // Sequential: 31 // 32 // ./area -problem cube -petscspace_degree 3 -dm_refine 2 33 // 34 // ./area -problem sphere -petscspace_degree 3 -dm_refine 2 35 // 36 // In parallel: 37 // 38 // mpiexec -n 4 ./area -probelm cube -petscspace_degree 3 -dm_refine 2 39 // 40 // mpiexec -n 4 ./area -problem sphere -petscspace_degree 3 -dm_refine 2 41 // 42 // The above example runs use 2 levels of refinement for the mesh. 43 // Use -dm_refine k, for k levels of uniform refinement. 44 // 45 //TESTARGS -ceed {ceed_resource} -test -petscspace_degree 3 46 47 /// @file 48 /// libCEED example using the mass operator to compute a cube or a cubed-sphere surface area using PETSc with DMPlex 49 static const char help[] = 50 "Compute surface area of a cube or a cubed-sphere using DMPlex in PETSc\n"; 51 52 #include <string.h> 53 #include <petscdmplex.h> 54 #include <ceed.h> 55 #include "qfunctions/area/areacube.h" 56 #include "qfunctions/area/areasphere.h" 57 58 #ifndef M_PI 59 # define M_PI 3.14159265358979323846 60 #endif 61 62 // Auxiliary function to define CEED restrictions from DMPlex data 63 static int CreateRestrictionPlex(Ceed ceed, CeedInterlaceMode imode, CeedInt P, 64 CeedInt ncomp, CeedElemRestriction *Erestrict, 65 DM dm) { 66 PetscInt ierr; 67 PetscInt c, cStart, cEnd, nelem, nnodes, *erestrict, eoffset; 68 PetscSection section; 69 Vec Uloc; 70 71 PetscFunctionBegin; 72 73 // Get Nelem 74 ierr = DMGetSection(dm, §ion); CHKERRQ(ierr); 75 ierr = DMPlexGetHeightStratum(dm, 0, &cStart,& cEnd); CHKERRQ(ierr); 76 nelem = cEnd - cStart; 77 78 // Get indices 79 ierr = PetscMalloc1(nelem*P*P, &erestrict); CHKERRQ(ierr); 80 for (c=cStart, eoffset = 0; c<cEnd; c++) { 81 PetscInt numindices, *indices, i; 82 ierr = DMPlexGetClosureIndices(dm, section, section, c, &numindices, 83 &indices, NULL); CHKERRQ(ierr); 84 for (i=0; i<numindices; i+=ncomp) { 85 for (PetscInt j=0; j<ncomp; j++) { 86 if (indices[i+j] != indices[i] + (PetscInt)(copysign(j, indices[i]))) 87 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, 88 "Cell %D closure indices not interlaced", c); 89 } 90 // NO BC on closed surfaces 91 PetscInt loc = indices[i]; 92 erestrict[eoffset++] = loc/ncomp; 93 } 94 ierr = DMPlexRestoreClosureIndices(dm, section, section, c, &numindices, 95 &indices, NULL); CHKERRQ(ierr); 96 } 97 98 // Setup CEED restriction 99 ierr = DMGetLocalVector(dm, &Uloc); CHKERRQ(ierr); 100 ierr = VecGetLocalSize(Uloc, &nnodes); CHKERRQ(ierr); 101 102 ierr = DMRestoreLocalVector(dm, &Uloc); CHKERRQ(ierr); 103 CeedElemRestrictionCreate(ceed, imode, nelem, P*P, nnodes/ncomp, ncomp, 104 CEED_MEM_HOST, CEED_COPY_VALUES, erestrict, 105 Erestrict); 106 ierr = PetscFree(erestrict); CHKERRQ(ierr); 107 108 PetscFunctionReturn(0); 109 } 110 111 // Utility function taken from petsc/src/dm/impls/plex/examples/tutorials/ex7.c 112 static PetscErrorCode ProjectToUnitSphere(DM dm) { 113 Vec coordinates; 114 PetscScalar *coords; 115 PetscInt Nv, v, dim, d; 116 PetscErrorCode ierr; 117 118 PetscFunctionBeginUser; 119 ierr = DMGetCoordinatesLocal(dm, &coordinates); CHKERRQ(ierr); 120 ierr = VecGetLocalSize(coordinates, &Nv); CHKERRQ(ierr); 121 ierr = VecGetBlockSize(coordinates, &dim); CHKERRQ(ierr); 122 Nv /= dim; 123 ierr = VecGetArray(coordinates, &coords); CHKERRQ(ierr); 124 for (v = 0; v < Nv; ++v) { 125 PetscReal r = 0.0; 126 127 for (d = 0; d < dim; ++d) r += PetscSqr(PetscRealPart(coords[v*dim+d])); 128 r = PetscSqrtReal(r); 129 for (d = 0; d < dim; ++d) coords[v*dim+d] /= r; 130 } 131 ierr = VecRestoreArray(coordinates, &coords); CHKERRQ(ierr); 132 PetscFunctionReturn(0); 133 } 134 135 int main(int argc, char **argv) { 136 PetscInt ierr; 137 MPI_Comm comm; 138 char filename[PETSC_MAX_PATH_LEN], 139 ceedresource[PETSC_MAX_PATH_LEN] = "/cpu/self"; 140 PetscInt lsize, gsize, xlsize, 141 qextra = 1, // default number of extra quadrature points 142 ncompx = 3, // number of components of 3D physical coordinates 143 ncompu = 1, // dimension of field to which apply mass operator 144 topodim = 2, // topological dimension of manifold 145 degree = 3; // default degree for finite element bases 146 PetscBool read_mesh = PETSC_FALSE, 147 test_mode = PETSC_FALSE; 148 PetscSpace sp; 149 PetscFE fe; 150 Vec X, Xloc, V, Vloc; 151 DM dm, dmcoord; 152 Ceed ceed; 153 CeedInt P, Q; 154 CeedOperator op_setupgeo, op_apply; 155 CeedQFunction qf_setupgeo, qf_apply; 156 CeedBasis basisx, basisu; 157 CeedInterlaceMode imodeceed = CEED_NONINTERLACED, 158 imodepetsc = CEED_INTERLACED; 159 CeedElemRestriction Erestrictx, Erestrictu, Erestrictxi, 160 Erestrictqdi; 161 PetscFunctionList geomfactorlist = NULL; 162 char problemtype[PETSC_MAX_PATH_LEN] = "sphere"; 163 164 ierr = PetscInitialize(&argc, &argv, NULL, help); 165 if (ierr) return ierr; 166 comm = PETSC_COMM_WORLD; 167 168 // Set up problem type command line option 169 ierr = PetscFunctionListAdd(&geomfactorlist, "cube", &SetupMassGeoCube); 170 CHKERRQ(ierr); 171 ierr = PetscFunctionListAdd(&geomfactorlist, "sphere", &SetupMassGeoSphere); 172 CHKERRQ(ierr); 173 174 // Read command line options 175 ierr = PetscOptionsBegin(comm, NULL, "CEED surface area problem with PETSc", 176 NULL); 177 CHKERRQ(ierr); 178 ierr = PetscOptionsFList("-problem", "Problem to solve", NULL, geomfactorlist, 179 problemtype, problemtype, sizeof problemtype, NULL); 180 CHKERRQ(ierr); 181 ierr = PetscOptionsInt("-qextra", "Number of extra quadrature points", 182 NULL, qextra, &qextra, NULL); CHKERRQ(ierr); 183 ierr = PetscOptionsString("-ceed", "CEED resource specifier", 184 NULL, ceedresource, ceedresource, 185 sizeof(ceedresource), NULL); CHKERRQ(ierr); 186 ierr = PetscOptionsBool("-test", 187 "Testing mode (do not print unless error is large)", 188 NULL, test_mode, &test_mode, NULL); CHKERRQ(ierr); 189 ierr = PetscOptionsString("-mesh", "Read mesh from file", NULL, 190 filename, filename, sizeof(filename), &read_mesh); 191 CHKERRQ(ierr); 192 ierr = PetscOptionsEnd(); CHKERRQ(ierr); 193 194 // Setup function pointer for geometric factors 195 int (*geomfp)(void *ctx, const CeedInt Q, const CeedScalar *const *in, 196 CeedScalar *const *out); 197 ierr = PetscFunctionListFind(geomfactorlist, problemtype, 198 (void(* *)(void))&geomfp); CHKERRQ(ierr); 199 const char *str; 200 if (geomfp == SetupMassGeoCube) 201 str = SetupMassGeoCube_loc; 202 else if (geomfp == SetupMassGeoSphere) 203 str = SetupMassGeoSphere_loc; 204 else 205 return CeedError(ceed, 1, "Function not found in the list"); 206 207 // Setup DM 208 if (read_mesh) { 209 ierr = DMPlexCreateFromFile(PETSC_COMM_WORLD, filename, PETSC_TRUE, &dm); 210 CHKERRQ(ierr); 211 } else { 212 // Create the mesh as a 0-refined sphere. This will create a cubic surface, not a box. 213 PetscBool simplex = PETSC_FALSE; 214 ierr = DMPlexCreateSphereMesh(PETSC_COMM_WORLD, topodim, simplex, &dm); 215 CHKERRQ(ierr); 216 // Set the object name 217 ierr = PetscObjectSetName((PetscObject)dm, problemtype); CHKERRQ(ierr); 218 // Distribute mesh over processes 219 { 220 DM dmDist = NULL; 221 PetscPartitioner part; 222 223 ierr = DMPlexGetPartitioner(dm, &part); CHKERRQ(ierr); 224 ierr = PetscPartitionerSetFromOptions(part); CHKERRQ(ierr); 225 ierr = DMPlexDistribute(dm, 0, NULL, &dmDist); CHKERRQ(ierr); 226 if (dmDist) { 227 ierr = DMDestroy(&dm); CHKERRQ(ierr); 228 dm = dmDist; 229 } 230 } 231 // Refine DMPlex with uniform refinement using runtime option -dm_refine 232 ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE); CHKERRQ(ierr); 233 ierr = DMSetFromOptions(dm); CHKERRQ(ierr); 234 if (!strcmp(problemtype, "sphere")) 235 ierr = ProjectToUnitSphere(dm); CHKERRQ(ierr); 236 // View DMPlex via runtime option 237 ierr = DMViewFromOptions(dm, NULL, "-dm_view"); CHKERRQ(ierr); 238 } 239 240 // Create FE 241 ierr = PetscFECreateDefault(PETSC_COMM_SELF, topodim, ncompu, PETSC_FALSE, NULL, 242 PETSC_DETERMINE, &fe); 243 CHKERRQ(ierr); 244 ierr = DMSetFromOptions(dm); CHKERRQ(ierr); 245 ierr = DMAddField(dm, NULL, (PetscObject)fe); CHKERRQ(ierr); 246 ierr = DMCreateDS(dm); CHKERRQ(ierr); 247 ierr = DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL); 248 CHKERRQ(ierr); 249 250 // Get basis space degree 251 ierr = PetscFEGetBasisSpace(fe, &sp); CHKERRQ(ierr); 252 ierr = PetscSpaceGetDegree(sp, °ree, NULL); CHKERRQ(ierr); 253 ierr = PetscFEDestroy(&fe); CHKERRQ(ierr); 254 if (degree < 1) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, 255 "-petscspace_degree %D must be at least 1", degree); 256 257 // Create vectors 258 ierr = DMCreateGlobalVector(dm, &X); CHKERRQ(ierr); 259 ierr = VecGetLocalSize(X, &lsize); CHKERRQ(ierr); 260 ierr = VecGetSize(X, &gsize); CHKERRQ(ierr); 261 ierr = DMCreateLocalVector(dm, &Xloc); CHKERRQ(ierr); 262 ierr = VecGetSize(Xloc, &xlsize); CHKERRQ(ierr); 263 ierr = VecDuplicate(X, &V); CHKERRQ(ierr); 264 ierr = VecDuplicate(Xloc, &Vloc); CHKERRQ(ierr); 265 266 // Set up libCEED 267 CeedInit(ceedresource, &ceed); 268 269 // Print summary 270 P = degree + 1; 271 Q = P + qextra; 272 const char *usedresource; 273 CeedGetResource(ceed, &usedresource); 274 if (!test_mode) { 275 ierr = PetscPrintf(comm, 276 "\n-- libCEED + PETSc Surface Area of a Manifold --\n" 277 " libCEED:\n" 278 " libCEED Backend : %s\n" 279 " Mesh:\n" 280 " Number of 1D Basis Nodes (p) : %d\n" 281 " Number of 1D Quadrature Points (q) : %d\n" 282 " Global nodes : %D\n" 283 " DoF per node : %D\n", 284 usedresource, P, Q, gsize/ncompu, ncompu); 285 CHKERRQ(ierr); 286 } 287 288 // Setup libCEED's objects: 289 // Create bases 290 CeedBasisCreateTensorH1Lagrange(ceed, topodim, ncompu, P, Q, 291 CEED_GAUSS, &basisu); 292 CeedBasisCreateTensorH1Lagrange(ceed, topodim, ncompx, 2, Q, 293 CEED_GAUSS, &basisx); 294 295 // CEED restrictions 296 ierr = DMGetCoordinateDM(dm, &dmcoord); CHKERRQ(ierr); 297 ierr = DMPlexSetClosurePermutationTensor(dmcoord, PETSC_DETERMINE, NULL); 298 CHKERRQ(ierr); 299 300 CreateRestrictionPlex(ceed, imodepetsc, 2, ncompx, &Erestrictx, dmcoord); 301 CHKERRQ(ierr); 302 CreateRestrictionPlex(ceed, imodepetsc, P, ncompu, &Erestrictu, dm); 303 CHKERRQ(ierr); 304 305 CeedInt cStart, cEnd; 306 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd); CHKERRQ(ierr); 307 const CeedInt nelem = cEnd - cStart; 308 309 // CEED identity restrictions 310 const CeedInt qdatasize = 1; 311 CeedElemRestrictionCreateIdentity(ceed, imodeceed, nelem, Q*Q, nelem*Q*Q, 312 qdatasize, &Erestrictqdi); 313 CeedElemRestrictionCreateIdentity(ceed, imodeceed, nelem, Q*Q, nelem*Q*Q, 1, 314 &Erestrictxi); 315 316 // Element coordinates 317 Vec coords; 318 const PetscScalar *coordArray; 319 PetscSection section; 320 ierr = DMGetCoordinatesLocal(dm, &coords); CHKERRQ(ierr); 321 ierr = VecGetArrayRead(coords, &coordArray); CHKERRQ(ierr); 322 ierr = DMGetSection(dmcoord, §ion); CHKERRQ(ierr); 323 324 CeedVector xcoord; 325 CeedElemRestrictionCreateVector(Erestrictx, &xcoord, NULL); 326 CeedVectorSetArray(xcoord, CEED_MEM_HOST, CEED_COPY_VALUES, 327 (PetscScalar *)coordArray); 328 ierr = VecRestoreArrayRead(coords, &coordArray); 329 330 // Create the vectors that will be needed in setup and apply 331 CeedVector uceed, vceed, qdata; 332 CeedInt nqpts; 333 CeedBasisGetNumQuadraturePoints(basisu, &nqpts); 334 CeedVectorCreate(ceed, qdatasize*nelem*nqpts, &qdata); 335 CeedVectorCreate(ceed, xlsize, &uceed); 336 CeedVectorCreate(ceed, xlsize, &vceed); 337 338 // Create the Q-function that builds the operator for the geomteric factors 339 // (i.e., the quadrature data) 340 CeedQFunctionCreateInterior(ceed, 1, geomfp, str, &qf_setupgeo); 341 CeedQFunctionAddInput(qf_setupgeo, "x", ncompx, CEED_EVAL_INTERP); 342 CeedQFunctionAddInput(qf_setupgeo, "dx", ncompx*topodim, CEED_EVAL_GRAD); 343 CeedQFunctionAddInput(qf_setupgeo, "weight", 1, CEED_EVAL_WEIGHT); 344 CeedQFunctionAddOutput(qf_setupgeo, "qdata", qdatasize, CEED_EVAL_NONE); 345 346 // Set up the mass operator 347 CeedQFunctionCreateInterior(ceed, 1, Mass, Mass_loc, &qf_apply); 348 CeedQFunctionAddInput(qf_apply, "u", ncompu, CEED_EVAL_INTERP); 349 CeedQFunctionAddInput(qf_apply, "qdata", qdatasize, CEED_EVAL_NONE); 350 CeedQFunctionAddOutput(qf_apply, "v", ncompu, CEED_EVAL_INTERP); 351 352 // Create the operator that builds the quadrature data for the operator 353 CeedOperatorCreate(ceed, qf_setupgeo, NULL, NULL, &op_setupgeo); 354 CeedOperatorSetField(op_setupgeo, "x", Erestrictx, basisx, 355 CEED_VECTOR_ACTIVE); 356 CeedOperatorSetField(op_setupgeo, "dx", Erestrictx, basisx, 357 CEED_VECTOR_ACTIVE); 358 CeedOperatorSetField(op_setupgeo, "weight", Erestrictxi, basisx, 359 CEED_VECTOR_NONE); 360 CeedOperatorSetField(op_setupgeo, "qdata", Erestrictqdi, 361 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 362 363 // Create the mass operator 364 CeedOperatorCreate(ceed, qf_apply, NULL, NULL, &op_apply); 365 CeedOperatorSetField(op_apply, "u", Erestrictu, basisu, CEED_VECTOR_ACTIVE); 366 CeedOperatorSetField(op_apply, "qdata", Erestrictqdi, CEED_BASIS_COLLOCATED, 367 qdata); 368 CeedOperatorSetField(op_apply, "v", Erestrictu, basisu, CEED_VECTOR_ACTIVE); 369 370 // Compute the quadrature data for the mass operator 371 CeedOperatorApply(op_setupgeo, xcoord, qdata, CEED_REQUEST_IMMEDIATE); 372 373 PetscScalar *v; 374 ierr = VecZeroEntries(Vloc); CHKERRQ(ierr); 375 ierr = VecGetArray(Vloc, &v); 376 CeedVectorSetArray(vceed, CEED_MEM_HOST, CEED_USE_POINTER, v); 377 378 // Compute the mesh volume using the mass operator: vol = 1^T \cdot M \cdot 1 379 if (!test_mode) { 380 ierr = PetscPrintf(comm, 381 "Computing the mesh area using the formula: area = 1^T M 1\n"); 382 CHKERRQ(ierr); 383 } 384 385 // Initialize u and v with ones 386 CeedVectorSetValue(uceed, 1.0); 387 388 // Apply the mass operator: 'u' -> 'v' 389 CeedOperatorApply(op_apply, uceed, vceed, CEED_REQUEST_IMMEDIATE); 390 CeedVectorSyncArray(vceed, CEED_MEM_HOST); 391 392 // Gather output vector 393 ierr = VecRestoreArray(Vloc, &v); CHKERRQ(ierr); 394 ierr = VecZeroEntries(V); CHKERRQ(ierr); 395 ierr = DMLocalToGlobalBegin(dm, Vloc, ADD_VALUES, V); CHKERRQ(ierr); 396 ierr = DMLocalToGlobalEnd(dm, Vloc, ADD_VALUES, V); CHKERRQ(ierr); 397 398 // Compute and print the sum of the entries of 'v' giving the mesh surface area 399 PetscScalar area; 400 ierr = VecSum(V, &area); CHKERRQ(ierr); 401 402 // Compute the exact surface area and print the result 403 CeedScalar exact_surfarea = 4 * M_PI; 404 if (!strcmp(problemtype, "cube")) { 405 PetscScalar l = 1.0/PetscSqrtReal(3.0); // half edge of the cube 406 exact_surfarea = 6 * (2*l) * (2*l); 407 } 408 409 if (!test_mode) { 410 ierr = PetscPrintf(comm, "Exact mesh surface area : % .14g\n", 411 exact_surfarea); CHKERRQ(ierr); 412 ierr = PetscPrintf(comm, "Computed mesh surface area : % .14g\n", area); 413 CHKERRQ(ierr); 414 ierr = PetscPrintf(comm, "Area error : % .14g\n", 415 fabs(area - exact_surfarea)); CHKERRQ(ierr); 416 } 417 418 // PETSc cleanup 419 ierr = DMDestroy(&dm); CHKERRQ(ierr); 420 ierr = VecDestroy(&X); CHKERRQ(ierr); 421 ierr = VecDestroy(&Xloc); CHKERRQ(ierr); 422 ierr = VecDestroy(&V); CHKERRQ(ierr); 423 ierr = VecDestroy(&Vloc); CHKERRQ(ierr); 424 425 // libCEED cleanup 426 CeedQFunctionDestroy(&qf_setupgeo); 427 CeedOperatorDestroy(&op_setupgeo); 428 CeedVectorDestroy(&xcoord); 429 CeedVectorDestroy(&uceed); 430 CeedVectorDestroy(&vceed); 431 CeedVectorDestroy(&qdata); 432 CeedBasisDestroy(&basisx); 433 CeedBasisDestroy(&basisu); 434 CeedElemRestrictionDestroy(&Erestrictu); 435 CeedElemRestrictionDestroy(&Erestrictx); 436 CeedElemRestrictionDestroy(&Erestrictxi); 437 CeedElemRestrictionDestroy(&Erestrictqdi); 438 CeedQFunctionDestroy(&qf_apply); 439 CeedOperatorDestroy(&op_apply); 440 CeedDestroy(&ceed); 441 return PetscFinalize(); 442 } 443