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