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