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