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