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