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