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 21 // 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 -petscspace_degree 3 33 // 34 // In parallel: 35 // 36 // mpiexec -n 4 area -petscspace_degree 3 37 // 38 //TESTARGS -ceed {ceed_resource} -test -petscspace_degree 3 39 40 /// @file 41 /// libCEED example using the mass operator to compute surface area using PETSc with DMPlex 42 static const char help[] = 43 "Compute surface area of a cube using DMPlex in PETSc\n"; 44 45 #include <string.h> 46 #include <petscdmplex.h> 47 #include <ceed.h> 48 #include "qfunctions/area/area.h" 49 50 // Auxiliary function to define CEED restrictions from DMPlex data 51 static int CreateRestrictionPlex(Ceed ceed, CeedInt P, CeedInt ncomp, 52 CeedElemRestriction *Erestrict, DM dm) { 53 PetscInt ierr; 54 PetscInt c, cStart, cEnd, nelem, nnodes, *erestrict, eoffset; 55 PetscSection section; 56 Vec Uloc; 57 58 PetscFunctionBegin; 59 60 // Get Nelem 61 ierr = DMGetSection(dm, §ion); CHKERRQ(ierr); 62 ierr = DMPlexGetHeightStratum(dm, 0, &cStart,& cEnd); CHKERRQ(ierr); 63 nelem = cEnd - cStart; 64 65 // Get indices 66 ierr = PetscMalloc1(nelem*P*P, &erestrict); CHKERRQ(ierr); 67 for (c=cStart, eoffset = 0; c<cEnd; c++) { 68 PetscInt numindices, *indices, i; 69 ierr = DMPlexGetClosureIndices(dm, section, section, c, &numindices, 70 &indices, NULL); CHKERRQ(ierr); 71 for (i=0; i<numindices; i+=ncomp) { 72 for (PetscInt j=0; j<ncomp; j++) { 73 if (indices[i+j] != indices[i] + (PetscInt)(copysign(j, indices[i]))) 74 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, 75 "Cell %D closure indices not interlaced", c); 76 } 77 // NO BC on closed surfaces 78 PetscInt loc = indices[i]; 79 erestrict[eoffset++] = loc/ncomp; 80 } 81 ierr = DMPlexRestoreClosureIndices(dm, section, section, c, &numindices, 82 &indices, NULL); CHKERRQ(ierr); 83 } 84 85 // Setup CEED restriction 86 ierr = DMGetLocalVector(dm, &Uloc); CHKERRQ(ierr); 87 ierr = VecGetLocalSize(Uloc, &nnodes); CHKERRQ(ierr); 88 89 ierr = DMRestoreLocalVector(dm, &Uloc); CHKERRQ(ierr); 90 CeedElemRestrictionCreate(ceed, nelem, P*P, nnodes/ncomp, ncomp, 91 CEED_MEM_HOST, CEED_COPY_VALUES, erestrict, 92 Erestrict); 93 ierr = PetscFree(erestrict); CHKERRQ(ierr); 94 95 PetscFunctionReturn(0); 96 } 97 98 int main(int argc, char **argv) { 99 PetscInt ierr; 100 MPI_Comm comm; 101 char filename[PETSC_MAX_PATH_LEN], 102 ceedresource[PETSC_MAX_PATH_LEN] = "/cpu/self"; 103 PetscInt lsize, gsize, xlsize, 104 qextra = 1, // default number of extra quadrature points 105 ncompx = 3, // number of components of 3D physical coordinates 106 ncompu = 1, // dimension of field to which apply mass operator 107 topodim = 2, // topological dimension of manifold 108 degree = 3; // default degree for finite element bases 109 PetscBool read_mesh = PETSC_FALSE, 110 test_mode = PETSC_FALSE; 111 PetscSpace sp; 112 PetscFE fe; 113 Vec X, Xloc, V, Vloc; 114 DM dm, dmcoord; 115 Ceed ceed; 116 CeedInt P, Q; 117 118 ierr = PetscInitialize(&argc, &argv, NULL, help); 119 if (ierr) return ierr; 120 comm = PETSC_COMM_WORLD; 121 122 // Read CL options 123 ierr = PetscOptionsBegin(comm, NULL, "CEED surface area problem with PETSc", 124 NULL); 125 CHKERRQ(ierr); 126 ierr = PetscOptionsInt("-qextra", "Number of extra quadrature points", 127 NULL, qextra, &qextra, NULL); CHKERRQ(ierr); 128 ierr = PetscOptionsString("-ceed", "CEED resource specifier", 129 NULL, ceedresource, ceedresource, 130 sizeof(ceedresource), NULL); CHKERRQ(ierr); 131 ierr = PetscOptionsBool("-test", 132 "Testing mode (do not print unless error is large)", 133 NULL, test_mode, &test_mode, NULL); CHKERRQ(ierr); 134 ierr = PetscOptionsString("-mesh", "Read mesh from file", NULL, 135 filename, filename, sizeof(filename), &read_mesh); 136 CHKERRQ(ierr); 137 ierr = PetscOptionsEnd(); CHKERRQ(ierr); 138 139 // Setup DM 140 PetscScalar l = 1.0/PetscSqrtReal(3.0); // half edge of the cube 141 if (read_mesh) { 142 ierr = DMPlexCreateFromFile(PETSC_COMM_WORLD, filename, PETSC_TRUE, &dm); 143 CHKERRQ(ierr); 144 } else { 145 // Create the mesh as a 0-refined sphere. This will create a cubic surface, not a box. 146 PetscBool simplex = PETSC_FALSE; 147 ierr = DMPlexCreateSphereMesh(PETSC_COMM_WORLD, topodim, simplex, &dm); 148 CHKERRQ(ierr); 149 // Set the object name 150 ierr = PetscObjectSetName((PetscObject) dm, "Cube"); CHKERRQ(ierr); 151 // Distribute mesh over processes 152 { 153 DM dmDist = NULL; 154 PetscPartitioner part; 155 156 ierr = DMPlexGetPartitioner(dm, &part); CHKERRQ(ierr); 157 ierr = PetscPartitionerSetFromOptions(part); CHKERRQ(ierr); 158 ierr = DMPlexDistribute(dm, 0, NULL, &dmDist); CHKERRQ(ierr); 159 if (dmDist) { 160 ierr = DMDestroy(&dm); CHKERRQ(ierr); 161 dm = dmDist; 162 } 163 } 164 // View DMPlex via runtime option 165 ierr = DMViewFromOptions(dm, NULL, "-dm_view"); CHKERRQ(ierr); 166 } 167 168 // Create FE 169 ierr = PetscFECreateDefault(PETSC_COMM_SELF, topodim, ncompu, PETSC_FALSE, NULL, 170 PETSC_DETERMINE, &fe); 171 CHKERRQ(ierr); 172 ierr = DMSetFromOptions(dm); CHKERRQ(ierr); 173 ierr = DMAddField(dm, NULL, (PetscObject)fe); CHKERRQ(ierr); 174 ierr = DMCreateDS(dm); CHKERRQ(ierr); 175 ierr = DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL); 176 CHKERRQ(ierr); 177 178 // Get basis space degree 179 ierr = PetscFEGetBasisSpace(fe, &sp); CHKERRQ(ierr); 180 ierr = PetscSpaceGetDegree(sp, °ree, NULL); CHKERRQ(ierr); 181 ierr = PetscFEDestroy(&fe); CHKERRQ(ierr); 182 if (degree < 1) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, 183 "-petscspace_degree %D must be at least 1", degree); 184 185 // Create vectors 186 ierr = DMCreateGlobalVector(dm, &X); CHKERRQ(ierr); 187 ierr = VecGetLocalSize(X, &lsize); CHKERRQ(ierr); 188 ierr = VecGetSize(X, &gsize); CHKERRQ(ierr); 189 ierr = DMCreateLocalVector(dm, &Xloc); CHKERRQ(ierr); 190 ierr = VecGetSize(Xloc, &xlsize); CHKERRQ(ierr); 191 ierr = VecDuplicate(X, &V); CHKERRQ(ierr); 192 ierr = VecDuplicate(Xloc, &Vloc); CHKERRQ(ierr); 193 194 // Set up libCEED 195 CeedInit(ceedresource, &ceed); 196 197 // Print summary 198 P = degree + 1; 199 Q = P + qextra; 200 const char *usedresource; 201 CeedGetResource(ceed, &usedresource); 202 if (!test_mode) { 203 ierr = PetscPrintf(comm, 204 "\n-- libCEED + PETSc Surface Area problem --\n" 205 " libCEED:\n" 206 " libCEED Backend : %s\n" 207 " Mesh:\n" 208 " Number of 1D Basis Nodes (p) : %d\n" 209 " Number of 1D Quadrature Points (q) : %d\n" 210 " Global nodes : %D\n", 211 usedresource, P, Q, gsize/ncompu); 212 CHKERRQ(ierr); 213 } 214 215 // Setup libCEED's objects 216 // Create CEED operators and API objects they need 217 CeedOperator op_setupgeo, op_apply; 218 CeedQFunction qf_setupgeo, qf_apply; 219 CeedBasis basisx, basisu; 220 CeedElemRestriction Erestrictx, Erestrictu, Erestrictxi, 221 Erestrictqdi; 222 223 // Create bases 224 CeedBasisCreateTensorH1Lagrange(ceed, topodim, ncompu, P, Q, 225 CEED_GAUSS, &basisu); 226 CeedBasisCreateTensorH1Lagrange(ceed, topodim, ncompx, 2, Q, 227 CEED_GAUSS, &basisx); 228 229 // CEED restrictions 230 ierr = DMGetCoordinateDM(dm, &dmcoord); CHKERRQ(ierr); 231 ierr = DMPlexSetClosurePermutationTensor(dmcoord, PETSC_DETERMINE, NULL); 232 CHKERRQ(ierr); 233 234 CreateRestrictionPlex(ceed, 2, ncompx, &Erestrictx, dmcoord); CHKERRQ(ierr); 235 CreateRestrictionPlex(ceed, P, ncompu, &Erestrictu, dm); CHKERRQ(ierr); 236 237 CeedInt cStart, cEnd; 238 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd); CHKERRQ(ierr); 239 const CeedInt nelem = cEnd - cStart; 240 241 // CEED identity restrictions 242 const CeedInt qdatasize = 1; 243 CeedElemRestrictionCreateIdentity(ceed, nelem, Q*Q, nelem*Q*Q, 244 qdatasize, &Erestrictqdi); 245 CeedElemRestrictionCreateIdentity(ceed, nelem, Q*Q, nelem*Q*Q, 1, 246 &Erestrictxi); 247 248 // Element coordinates 249 Vec coords; 250 const PetscScalar *coordArray; 251 PetscSection section; 252 ierr = DMGetCoordinatesLocal(dm, &coords); CHKERRQ(ierr); 253 ierr = VecGetArrayRead(coords, &coordArray); CHKERRQ(ierr); 254 ierr = DMGetSection(dmcoord, §ion); CHKERRQ(ierr); 255 256 CeedVector xcoord; 257 CeedElemRestrictionCreateVector(Erestrictx, &xcoord, NULL); 258 CeedVectorSetArray(xcoord, CEED_MEM_HOST, CEED_COPY_VALUES, 259 (PetscScalar *)coordArray); 260 ierr = VecRestoreArrayRead(coords, &coordArray); 261 262 // Create the vectors that will be needed in setup and apply 263 CeedVector uceed, vceed, qdata; 264 CeedInt nqpts; 265 CeedBasisGetNumQuadraturePoints(basisu, &nqpts); 266 CeedVectorCreate(ceed, qdatasize*nelem*nqpts, &qdata); 267 CeedVectorCreate(ceed, xlsize, &uceed); 268 CeedVectorCreate(ceed, xlsize, &vceed); 269 270 /* Create the Q-function that builds the operator for the geomteric factors 271 (i.e., the quadrature data) */ 272 CeedQFunctionCreateInterior(ceed, 1, SetupMassGeo, 273 SetupMassGeo_loc, &qf_setupgeo); 274 CeedQFunctionAddInput(qf_setupgeo, "dx", ncompx*topodim, CEED_EVAL_GRAD); 275 CeedQFunctionAddInput(qf_setupgeo, "weight", 1, CEED_EVAL_WEIGHT); 276 CeedQFunctionAddOutput(qf_setupgeo, "qdata", qdatasize, CEED_EVAL_NONE); 277 278 // Set up the mass operator 279 CeedQFunctionCreateInterior(ceed, 1, Mass, Mass_loc, &qf_apply); 280 CeedQFunctionAddInput(qf_apply, "u", ncompu, CEED_EVAL_INTERP); 281 CeedQFunctionAddInput(qf_apply, "qdata", qdatasize, CEED_EVAL_NONE); 282 CeedQFunctionAddOutput(qf_apply, "v", ncompu, CEED_EVAL_INTERP); 283 284 // Create the operator that builds the quadrature data for the operator 285 CeedOperatorCreate(ceed, qf_setupgeo, CEED_QFUNCTION_NONE, 286 CEED_QFUNCTION_NONE, &op_setupgeo); 287 CeedOperatorSetField(op_setupgeo, "dx", Erestrictx, CEED_TRANSPOSE, 288 basisx, CEED_VECTOR_ACTIVE); 289 CeedOperatorSetField(op_setupgeo, "weight", Erestrictxi, CEED_NOTRANSPOSE, 290 basisx, CEED_VECTOR_NONE); 291 CeedOperatorSetField(op_setupgeo, "qdata", Erestrictqdi, CEED_NOTRANSPOSE, 292 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 293 294 // Create the mass operator 295 CeedOperatorCreate(ceed, qf_apply, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 296 &op_apply); 297 CeedOperatorSetField(op_apply, "u", Erestrictu, CEED_TRANSPOSE, 298 basisu, CEED_VECTOR_ACTIVE); 299 CeedOperatorSetField(op_apply, "qdata", Erestrictqdi, CEED_NOTRANSPOSE, 300 CEED_BASIS_COLLOCATED, qdata); 301 CeedOperatorSetField(op_apply, "v", Erestrictu, CEED_TRANSPOSE, 302 basisu, CEED_VECTOR_ACTIVE); 303 304 // Compute the quadrature data for the mass operator 305 CeedOperatorApply(op_setupgeo, xcoord, qdata, CEED_REQUEST_IMMEDIATE); 306 307 PetscScalar *v; 308 ierr = VecZeroEntries(Vloc); CHKERRQ(ierr); 309 ierr = VecGetArray(Vloc, &v); 310 CeedVectorSetArray(vceed, CEED_MEM_HOST, CEED_USE_POINTER, v); 311 312 // Compute the mesh volume using the mass operator: vol = 1^T \cdot M \cdot 1 313 if (!test_mode) { 314 ierr = PetscPrintf(comm, 315 "Computing the mesh volume using the formula: vol = 1^T M 1\n"); 316 CHKERRQ(ierr); 317 } 318 319 // Initialize u and v with ones 320 CeedVectorSetValue(uceed, 1.0); 321 CeedVectorSetValue(vceed, 1.0); 322 323 // Apply the mass operator: 'u' -> 'v' 324 CeedOperatorApply(op_apply, uceed, vceed, CEED_REQUEST_IMMEDIATE); 325 CeedVectorSyncArray(vceed, CEED_MEM_HOST); 326 327 // Gather output vector 328 ierr = VecRestoreArray(Vloc, &v); CHKERRQ(ierr); 329 ierr = VecZeroEntries(V); CHKERRQ(ierr); 330 ierr = DMLocalToGlobalBegin(dm, Vloc, ADD_VALUES, V); CHKERRQ(ierr); 331 ierr = DMLocalToGlobalEnd(dm, Vloc, ADD_VALUES, V); CHKERRQ(ierr); 332 333 // Compute and print the sum of the entries of 'v' giving the mesh surface area 334 PetscScalar area; 335 ierr = VecSum(V, &area); CHKERRQ(ierr); 336 337 // Compute the exact surface area and print the result 338 CeedScalar exact_surfarea = 6 * (2*l) * (2*l); 339 if (!test_mode) { 340 ierr = PetscPrintf(comm, "Exact mesh surface area : % .14g\n", 341 exact_surfarea); CHKERRQ(ierr); 342 ierr = PetscPrintf(comm, "Computed mesh surface area : % .14g\n", area); 343 CHKERRQ(ierr); 344 ierr = PetscPrintf(comm, "Area error : % .14g\n", 345 fabs(area - exact_surfarea)); CHKERRQ(ierr); 346 } 347 348 // PETSc cleanup 349 ierr = DMDestroy(&dm); CHKERRQ(ierr); 350 ierr = VecDestroy(&X); CHKERRQ(ierr); 351 ierr = VecDestroy(&Xloc); CHKERRQ(ierr); 352 ierr = VecDestroy(&V); CHKERRQ(ierr); 353 ierr = VecDestroy(&Vloc); CHKERRQ(ierr); 354 355 // libCEED cleanup 356 CeedQFunctionDestroy(&qf_setupgeo); 357 CeedOperatorDestroy(&op_setupgeo); 358 CeedVectorDestroy(&xcoord); 359 CeedVectorDestroy(&uceed); 360 CeedVectorDestroy(&vceed); 361 CeedVectorDestroy(&qdata); 362 CeedBasisDestroy(&basisx); 363 CeedBasisDestroy(&basisu); 364 CeedElemRestrictionDestroy(&Erestrictu); 365 CeedElemRestrictionDestroy(&Erestrictx); 366 CeedElemRestrictionDestroy(&Erestrictxi); 367 CeedElemRestrictionDestroy(&Erestrictqdi); 368 CeedQFunctionDestroy(&qf_apply); 369 CeedOperatorDestroy(&op_apply); 370 CeedDestroy(&ceed); 371 return PetscFinalize(); 372 } 373