1 // Copyright (c) 2017-2018, Lawrence Livermore National Security, LLC. 2 // Produced at the Lawrence Livermore National Laboratory. LLNL-CODE-734707. 3 // All Rights 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 Example 2 18 // 19 // This example illustrates a simple usage of libCEED to compute the surface 20 // area of a 3D body using matrix-free application of a diffusion operator. 21 // Arbitrary mesh and solution orders in 1D, 2D and 3D are supported from the 22 // same code. 23 // 24 // The example has no dependencies, and is designed to be self-contained. For 25 // additional examples that use external discretization libraries (MFEM, PETSc, 26 // etc.) see the subdirectories in libceed/examples. 27 // 28 // All libCEED objects use a Ceed device object constructed based on a command 29 // line argument (-ceed). 30 // 31 // Build with: 32 // 33 // make ex2-surface [CEED_DIR=</path/to/libceed>] 34 // 35 // Sample runs: 36 // 37 // ./ex2-surface 38 // ./ex2-surface -ceed /cpu/self 39 // ./ex2-surface -ceed /gpu/occa 40 // ./ex2-surface -ceed /cpu/occa 41 // ./ex2-surface -ceed /omp/occa 42 // ./ex2-surface -ceed /ocl/occa 43 // ./ex2-surface -m ../../../mfem/data/fichera.mesh 44 // ./ex2-surface -m ../../../mfem/data/star.vtk -o 3 45 // ./ex2-surface -m ../../../mfem/data/inline-segment.mesh -o 8 46 // 47 // Next line is grep'd from tap.sh to set its arguments 48 // Test in 1D-3D 49 //TESTARGS -ceed {ceed_resource} -d 2 -t 50 //TESTARGS -ceed {ceed_resource} -d 3 -t 51 //TESTARGS -ceed {ceed_resource} -d 1 -t -g 52 //TESTARGS -ceed {ceed_resource} -d 3 -t -g 53 54 /// @file 55 /// libCEED example using diffusion operator to compute surface area 56 57 #include <ceed.h> 58 #include <stdlib.h> 59 #include <math.h> 60 #include <string.h> 61 62 #include "ex2-surface.h" 63 64 // Auxiliary functions. 65 int GetCartesianMeshSize(int dim, int order, int prob_size, int nxyz[3]); 66 int BuildCartesianRestriction(Ceed ceed, int dim, int nxyz[3], int order, 67 int ncomp, CeedInt *size, CeedInt num_qpts, 68 CeedElemRestriction *restr, 69 CeedElemRestriction *restr_i); 70 int SetCartesianMeshCoords(int dim, int nxyz[3], int mesh_order, 71 CeedVector mesh_coords); 72 CeedScalar TransformMeshCoords(int dim, int mesh_size, CeedVector mesh_coords); 73 74 75 int main(int argc, const char *argv[]) { 76 const char *ceed_spec = "/cpu/self"; 77 int dim = 3; // dimension of the mesh 78 int ncompx = 3; // number of x components 79 int mesh_order = 4; // polynomial degree for the mesh 80 int sol_order = 4; // polynomial degree for the solution 81 int num_qpts = sol_order+2; // number of 1D quadrature points 82 int prob_size = -1; // approximate problem size 83 int help = 0, test = 0, gallery = 0; 84 85 // Process command line arguments. 86 for (int ia = 1; ia < argc; ia++) { 87 int next_arg = ((ia+1) < argc), parse_error = 0; 88 if (!strcmp(argv[ia],"-h")) { 89 help = 1; 90 } else if (!strcmp(argv[ia],"-c") || !strcmp(argv[ia],"-ceed")) { 91 parse_error = next_arg ? ceed_spec = argv[++ia], 0 : 1; 92 } else if (!strcmp(argv[ia],"-d")) { 93 parse_error = next_arg ? dim = atoi(argv[++ia]), 0 : 1; 94 ncompx = dim; 95 } else if (!strcmp(argv[ia],"-m")) { 96 parse_error = next_arg ? mesh_order = atoi(argv[++ia]), 0 : 1; 97 } else if (!strcmp(argv[ia],"-o")) { 98 parse_error = next_arg ? sol_order = atoi(argv[++ia]), 0 : 1; 99 } else if (!strcmp(argv[ia],"-q")) { 100 parse_error = next_arg ? num_qpts = atoi(argv[++ia]), 0 : 1; 101 } else if (!strcmp(argv[ia],"-s")) { 102 parse_error = next_arg ? prob_size = atoi(argv[++ia]), 0 : 1; 103 } else if (!strcmp(argv[ia],"-t")) { 104 test = 1; 105 } else if (!strcmp(argv[ia],"-g")) { 106 gallery = 1; 107 } 108 if (parse_error) { 109 printf("Error parsing command line options.\n"); 110 return 1; 111 } 112 } 113 if (prob_size < 0) prob_size = test ? 16*16*dim*dim : 256*1024; 114 115 // Set mesh_order = sol_order. 116 mesh_order = fmax(mesh_order, sol_order); 117 sol_order = mesh_order; 118 119 // Print the values of all options: 120 if (!test || help) { 121 printf("Selected options: [command line option] : <current value>\n"); 122 printf(" Ceed specification [-c] : %s\n", ceed_spec); 123 printf(" Mesh dimension [-d] : %d\n", dim); 124 printf(" Mesh order [-m] : %d\n", mesh_order); 125 printf(" Solution order [-o] : %d\n", sol_order); 126 printf(" Num. 1D quadr. pts [-q] : %d\n", num_qpts); 127 printf(" Approx. # unknowns [-s] : %d\n", prob_size); 128 printf(" QFunction source [-g] : %s\n", gallery?"gallery":"header"); 129 if (help) { 130 printf("Test/quiet mode is %s\n", (test?"ON":"OFF (use -t to enable)")); 131 return 0; 132 } 133 printf("\n"); 134 } 135 136 // Select appropriate backend and logical device based on the <ceed-spec> 137 // command line argument. 138 Ceed ceed; 139 CeedInit(ceed_spec, &ceed); 140 141 // Construct the mesh and solution bases. 142 CeedBasis mesh_basis, sol_basis; 143 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompx, mesh_order+1, num_qpts, 144 CEED_GAUSS, &mesh_basis); 145 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, sol_order+1, num_qpts, 146 CEED_GAUSS, &sol_basis); 147 148 // Determine the mesh size based on the given approximate problem size. 149 int nxyz[3]; 150 GetCartesianMeshSize(dim, sol_order, prob_size, nxyz); 151 152 if (!test) { 153 printf("Mesh size: nx = %d", nxyz[0]); 154 if (dim > 1) { printf(", ny = %d", nxyz[1]); } 155 if (dim > 2) { printf(", nz = %d", nxyz[2]); } 156 printf("\n"); 157 } 158 159 // Build CeedElemRestriction objects describing the mesh and solution discrete 160 // representations. 161 CeedInt mesh_size, sol_size; 162 CeedElemRestriction mesh_restr, sol_restr, sol_restr_i, qdata_restr_i; 163 BuildCartesianRestriction(ceed, dim, nxyz, mesh_order, ncompx, &mesh_size, 164 num_qpts, &mesh_restr, NULL); 165 BuildCartesianRestriction(ceed, dim, nxyz, sol_order, dim*(dim+1)/2, 166 &sol_size, num_qpts, NULL, &qdata_restr_i); 167 BuildCartesianRestriction(ceed, dim, nxyz, sol_order, 1, &sol_size, 168 num_qpts, &sol_restr, &sol_restr_i); 169 if (!test) { 170 printf("Number of mesh nodes : %d\n", mesh_size/dim); 171 printf("Number of solution nodes : %d\n", sol_size); 172 } 173 174 // Create a CeedVector with the mesh coordinates. 175 CeedVector mesh_coords; 176 CeedVectorCreate(ceed, mesh_size, &mesh_coords); 177 SetCartesianMeshCoords(dim, nxyz, mesh_order, mesh_coords); 178 179 // Apply a transformation to the mesh. 180 CeedScalar exact_sa = TransformMeshCoords(dim, mesh_size, mesh_coords); 181 182 // Context data to be passed to the 'f_build_diff' Q-function. 183 CeedQFunctionContext build_ctx; 184 struct BuildContext build_ctx_data; 185 build_ctx_data.dim = build_ctx_data.space_dim = dim; 186 CeedQFunctionContextCreate(ceed, &build_ctx); 187 CeedQFunctionContextSetData(build_ctx, CEED_MEM_HOST, CEED_USE_POINTER, 188 sizeof(build_ctx_data), &build_ctx_data); 189 190 // Create the Q-function that builds the diffusion operator (i.e. computes its 191 // quadrature data) and set its context data. 192 CeedQFunction build_qfunc; 193 switch (gallery) { 194 case 0: 195 // This creates the QFunction directly. 196 CeedQFunctionCreateInterior(ceed, 1, f_build_diff, 197 f_build_diff_loc, &build_qfunc); 198 CeedQFunctionAddInput(build_qfunc, "dx", ncompx*dim, CEED_EVAL_GRAD); 199 CeedQFunctionAddInput(build_qfunc, "weights", 1, CEED_EVAL_WEIGHT); 200 CeedQFunctionAddOutput(build_qfunc, "qdata", dim*(dim+1)/2, CEED_EVAL_NONE); 201 CeedQFunctionSetContext(build_qfunc, build_ctx); 202 break; 203 case 1: { 204 // This creates the QFunction via the gallery. 205 char name[16] = ""; 206 snprintf(name, sizeof name, "Poisson%dDBuild", dim); 207 CeedQFunctionCreateInteriorByName(ceed, name, &build_qfunc); 208 break; 209 } 210 } 211 212 // Create the operator that builds the quadrature data for the diffusion 213 // operator. 214 CeedOperator build_oper; 215 CeedOperatorCreate(ceed, build_qfunc, CEED_QFUNCTION_NONE, 216 CEED_QFUNCTION_NONE, &build_oper); 217 CeedOperatorSetField(build_oper, "dx", mesh_restr, mesh_basis, 218 CEED_VECTOR_ACTIVE); 219 CeedOperatorSetField(build_oper, "weights", CEED_ELEMRESTRICTION_NONE, 220 mesh_basis, CEED_VECTOR_NONE); 221 CeedOperatorSetField(build_oper, "qdata", qdata_restr_i, 222 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 223 224 // Compute the quadrature data for the diffusion operator. 225 CeedVector qdata; 226 CeedInt elem_qpts = CeedIntPow(num_qpts, dim); 227 CeedInt num_elem = 1; 228 for (int d = 0; d < dim; d++) 229 num_elem *= nxyz[d]; 230 CeedVectorCreate(ceed, num_elem*elem_qpts*dim*(dim+1)/2, &qdata); 231 if (!test) { 232 printf("Computing the quadrature data for the diffusion operator ..."); 233 fflush(stdout); 234 } 235 CeedOperatorApply(build_oper, mesh_coords, qdata, 236 CEED_REQUEST_IMMEDIATE); 237 if (!test) { 238 printf(" done.\n"); 239 } 240 241 // Create the Q-function that defines the action of the diffusion operator. 242 CeedQFunction apply_qfunc; 243 switch (gallery) { 244 case 0: 245 // This creates the QFunction directly. 246 CeedQFunctionCreateInterior(ceed, 1, f_apply_diff, 247 f_apply_diff_loc, &apply_qfunc); 248 CeedQFunctionAddInput(apply_qfunc, "du", dim, CEED_EVAL_GRAD); 249 CeedQFunctionAddInput(apply_qfunc, "qdata", dim*(dim+1)/2, CEED_EVAL_NONE); 250 CeedQFunctionAddOutput(apply_qfunc, "dv", dim, CEED_EVAL_GRAD); 251 CeedQFunctionSetContext(apply_qfunc, build_ctx); 252 break; 253 case 1: { 254 // This creates the QFunction via the gallery. 255 char name[16] = ""; 256 snprintf(name, sizeof name, "Poisson%dDApply", dim); 257 CeedQFunctionCreateInteriorByName(ceed, name, &apply_qfunc); 258 break; 259 } 260 } 261 262 // Create the diffusion operator. 263 CeedOperator oper; 264 CeedOperatorCreate(ceed, apply_qfunc, CEED_QFUNCTION_NONE, 265 CEED_QFUNCTION_NONE, &oper); 266 CeedOperatorSetField(oper, "du", sol_restr, sol_basis, CEED_VECTOR_ACTIVE); 267 CeedOperatorSetField(oper, "qdata", qdata_restr_i, CEED_BASIS_COLLOCATED, 268 qdata); 269 CeedOperatorSetField(oper, "dv", sol_restr, sol_basis, CEED_VECTOR_ACTIVE); 270 271 // Compute the mesh surface area using the diff operator: 272 // sa = 1^T \cdot abs( K \cdot x). 273 if (!test) { 274 printf("Computing the mesh surface area using the formula: sa = 1^T.|K.x| ..."); 275 fflush(stdout); 276 } 277 278 // Create auxiliary solution-size vectors. 279 CeedVector u, v; 280 CeedVectorCreate(ceed, sol_size, &u); 281 CeedVectorCreate(ceed, sol_size, &v); 282 283 // Initialize 'u' with sum of coordinates, x+y+z. 284 CeedScalar *u_host; 285 const CeedScalar *x_host; 286 CeedVectorGetArray(u, CEED_MEM_HOST, &u_host); 287 CeedVectorGetArrayRead(mesh_coords, CEED_MEM_HOST, &x_host); 288 for (CeedInt i = 0; i < sol_size; i++) { 289 u_host[i] = 0; 290 for (CeedInt d = 0; d < dim; d++) 291 u_host[i] += x_host[i+d*sol_size]; 292 } 293 CeedVectorRestoreArray(u, &u_host); 294 CeedVectorRestoreArrayRead(mesh_coords, &x_host); 295 296 // Apply the diffusion operator: 'u' -> 'v'. 297 CeedOperatorApply(oper, u, v, CEED_REQUEST_IMMEDIATE); 298 299 // Compute and print the sum of the entries of 'v' giving the mesh surface area. 300 const CeedScalar *v_host; 301 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_host); 302 CeedScalar sa = 0.; 303 for (CeedInt i = 0; i < sol_size; i++) { 304 sa += fabs(v_host[i]); 305 } 306 CeedVectorRestoreArrayRead(v, &v_host); 307 if (!test) { 308 printf(" done.\n"); 309 printf("Exact mesh surface area : % .14g\n", exact_sa); 310 printf("Computed mesh surface area : % .14g\n", sa); 311 printf("Surface area error : % .14g\n", sa-exact_sa); 312 } else { 313 CeedScalar tol = (dim==1? 1E-12 : dim==2? 1E-1 : 1E-1); 314 if (fabs(sa-exact_sa)>tol) 315 printf("Surface area error : % .14g\n", sa-exact_sa); 316 } 317 318 // Free dynamically allocated memory. 319 CeedVectorDestroy(&u); 320 CeedVectorDestroy(&v); 321 CeedVectorDestroy(&qdata); 322 CeedVectorDestroy(&mesh_coords); 323 CeedOperatorDestroy(&oper); 324 CeedQFunctionDestroy(&apply_qfunc); 325 CeedQFunctionContextDestroy(&build_ctx); 326 CeedOperatorDestroy(&build_oper); 327 CeedQFunctionDestroy(&build_qfunc); 328 CeedElemRestrictionDestroy(&sol_restr); 329 CeedElemRestrictionDestroy(&mesh_restr); 330 CeedElemRestrictionDestroy(&sol_restr_i); 331 CeedElemRestrictionDestroy(&qdata_restr_i); 332 CeedBasisDestroy(&sol_basis); 333 CeedBasisDestroy(&mesh_basis); 334 CeedDestroy(&ceed); 335 return 0; 336 } 337 338 339 int GetCartesianMeshSize(int dim, int order, int prob_size, int nxyz[3]) { 340 // Use the approximate formula: 341 // prob_size ~ num_elem * order^dim 342 CeedInt num_elem = prob_size / CeedIntPow(order, dim); 343 CeedInt s = 0; // find s: num_elem/2 < 2^s <= num_elem 344 while (num_elem > 1) { 345 num_elem /= 2; 346 s++; 347 } 348 CeedInt r = s%dim; 349 for (int d = 0; d < dim; d++) { 350 int sd = s/dim; 351 if (r > 0) { sd++; r--; } 352 nxyz[d] = 1 << sd; 353 } 354 return 0; 355 } 356 357 int BuildCartesianRestriction(Ceed ceed, int dim, int nxyz[3], int order, 358 int ncomp, CeedInt *size, CeedInt num_qpts, 359 CeedElemRestriction *restr, 360 CeedElemRestriction *restr_i) { 361 CeedInt p = order, pp1 = p+1; 362 CeedInt nnodes = CeedIntPow(pp1, dim); // number of scal. nodes per element 363 CeedInt elem_qpts = CeedIntPow(num_qpts, dim); // number of qpts per element 364 CeedInt nd[3], num_elem = 1, scalar_size = 1; 365 for (int d = 0; d < dim; d++) { 366 num_elem *= nxyz[d]; 367 nd[d] = nxyz[d]*p + 1; 368 scalar_size *= nd[d]; 369 } 370 *size = scalar_size*ncomp; 371 // elem: 0 1 n-1 372 // |---*-...-*---|---*-...-*---|- ... -|--...--| 373 // nnodes: 0 1 p-1 p p+1 2*p n*p 374 CeedInt *el_nodes = malloc(sizeof(CeedInt)*num_elem*nnodes); 375 for (CeedInt e = 0; e < num_elem; e++) { 376 CeedInt exyz[3] = {1, 1, 1}, re = e; 377 for (int d = 0; d < dim; d++) { exyz[d] = re%nxyz[d]; re /= nxyz[d]; } 378 CeedInt *loc_el_nodes = el_nodes + e*nnodes; 379 for (int lnodes = 0; lnodes < nnodes; lnodes++) { 380 CeedInt gnodes = 0, gnodes_stride = 1, rnodes = lnodes; 381 for (int d = 0; d < dim; d++) { 382 gnodes += (exyz[d]*p + rnodes%pp1) * gnodes_stride; 383 gnodes_stride *= nd[d]; 384 rnodes /= pp1; 385 } 386 loc_el_nodes[lnodes] = gnodes; 387 } 388 } 389 if (restr) 390 CeedElemRestrictionCreate(ceed, num_elem, nnodes, ncomp, scalar_size, 391 ncomp*scalar_size, CEED_MEM_HOST, 392 CEED_COPY_VALUES, el_nodes, restr); 393 free(el_nodes); 394 395 if (restr_i) { 396 CeedElemRestrictionCreateStrided(ceed, num_elem, elem_qpts, 397 ncomp, ncomp*elem_qpts*num_elem, 398 CEED_STRIDES_BACKEND, restr_i); 399 } 400 401 return 0; 402 } 403 404 int SetCartesianMeshCoords(int dim, int nxyz[3], int mesh_order, 405 CeedVector mesh_coords) { 406 CeedInt p = mesh_order; 407 CeedInt nd[3], num_elem = 1, scalar_size = 1; 408 for (int d = 0; d < dim; d++) { 409 num_elem *= nxyz[d]; 410 nd[d] = nxyz[d]*p + 1; 411 scalar_size *= nd[d]; 412 } 413 CeedScalar *coords; 414 CeedVectorGetArray(mesh_coords, CEED_MEM_HOST, &coords); 415 CeedScalar *nodes = malloc(sizeof(CeedScalar)*(p+1)); 416 // The H1 basis uses Lobatto quadrature points as nodes. 417 CeedLobattoQuadrature(p+1, nodes, NULL); // nodes are in [-1,1] 418 for (CeedInt i = 0; i <= p; i++) { nodes[i] = 0.5+0.5*nodes[i]; } 419 for (CeedInt gsnodes = 0; gsnodes < scalar_size; gsnodes++) { 420 CeedInt rnodes = gsnodes; 421 for (int d = 0; d < dim; d++) { 422 CeedInt d1d = rnodes%nd[d]; 423 coords[gsnodes+scalar_size*d] = ((d1d/p)+nodes[d1d%p]) / nxyz[d]; 424 rnodes /= nd[d]; 425 } 426 } 427 free(nodes); 428 CeedVectorRestoreArray(mesh_coords, &coords); 429 return 0; 430 } 431 432 #ifndef M_PI 433 #define M_PI 3.14159265358979323846 434 #endif 435 436 CeedScalar TransformMeshCoords(int dim, int mesh_size, CeedVector mesh_coords) { 437 CeedScalar exact_sa = (dim==1? 2 : dim==2? 4 : 6); 438 CeedScalar *coords; 439 440 CeedVectorGetArray(mesh_coords, CEED_MEM_HOST, &coords); 441 for (CeedInt i = 0; i < mesh_size; i++) { 442 // map [0,1] to [0,1] varying the mesh density 443 coords[i] = 0.5+1./sqrt(3.)*sin((2./3.)*M_PI*(coords[i]-0.5)); 444 } 445 CeedVectorRestoreArray(mesh_coords, &coords); 446 447 return exact_sa; 448 } 449