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