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