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