1 // Copyright (c) 2017-2026, 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
main(int argc,const char * argv[])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, benchmark = 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], "-b")) {
85 parse_error = next_arg ? benchmark = 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 ? 16 * 16 * dim * dim : 256 * 1024;
98
99 // Set mesh_degree = sol_degree.
100 mesh_degree = fmax(mesh_degree, sol_degree);
101 sol_degree = mesh_degree;
102
103 // Print the values of all options:
104 if (!test || help) {
105 // LCOV_EXCL_START
106 printf("Selected options: [command line option] : <current value>\n");
107 printf(" Ceed specification [-c] : %s\n", ceed_spec);
108 printf(" Mesh dimension [-d] : %" CeedInt_FMT "\n", dim);
109 printf(" Mesh degree [-m] : %" CeedInt_FMT "\n", mesh_degree);
110 printf(" Solution degree [-p] : %" CeedInt_FMT "\n", sol_degree);
111 printf(" Num. 1D quadrature pts [-q] : %" CeedInt_FMT "\n", num_qpts);
112 printf(" Approx. # unknowns [-s] : %" CeedInt_FMT "\n", prob_size);
113 printf(" QFunction source [-g] : %s\n", gallery ? "gallery" : "header");
114 if (help) {
115 printf("Test/quiet mode is %s\n", (test ? "ON" : "OFF (use -t to enable)"));
116 return 0;
117 }
118 printf("\n");
119 // LCOV_EXCL_STOP
120 }
121
122 // Select appropriate backend and logical device based on the (-ceed) command line argument.
123 Ceed ceed;
124
125 CeedInit(ceed_spec, &ceed);
126
127 // Construct the mesh and solution bases.
128 CeedBasis mesh_basis, sol_basis;
129
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
136 GetCartesianMeshSize(dim, sol_degree, prob_size, num_xyz);
137
138 if (!test) {
139 // LCOV_EXCL_START
140 printf("Mesh size: nx = %" CeedInt_FMT, num_xyz[0]);
141 if (dim > 1) printf(", ny = %" CeedInt_FMT, num_xyz[1]);
142 if (dim > 2) printf(", nz = %" CeedInt_FMT, num_xyz[2]);
143 printf("\n");
144 // LCOV_EXCL_STOP
145 }
146
147 // Build CeedElemRestriction objects describing the mesh and solution discrete representations.
148 CeedInt mesh_size, sol_size;
149 CeedElemRestriction mesh_restriction, sol_restriction, q_data_restriction;
150
151 BuildCartesianRestriction(ceed, dim, num_xyz, mesh_degree, num_comp_x, &mesh_size, num_qpts, &mesh_restriction, NULL);
152 BuildCartesianRestriction(ceed, dim, num_xyz, sol_degree, dim * (dim + 1) / 2, &sol_size, num_qpts, NULL, &q_data_restriction);
153 BuildCartesianRestriction(ceed, dim, num_xyz, sol_degree, 1, &sol_size, num_qpts, &sol_restriction, NULL);
154 if (!test) {
155 // LCOV_EXCL_START
156 printf("Number of mesh nodes : %" CeedInt_FMT "\n", mesh_size / dim);
157 printf("Number of solution nodes : %" CeedInt_FMT "\n", sol_size);
158 // LCOV_EXCL_STOP
159 }
160
161 // Create a CeedVector with the mesh coordinates.
162 CeedVector mesh_coords;
163
164 CeedVectorCreate(ceed, mesh_size, &mesh_coords);
165 SetCartesianMeshCoords(dim, num_xyz, mesh_degree, mesh_coords);
166
167 // Apply a transformation to the mesh.
168 CeedScalar exact_surface_area = TransformMeshCoords(dim, mesh_size, mesh_coords);
169
170 // Context data to be passed to the 'build_diff' QFunction.
171 CeedQFunctionContext build_ctx;
172 struct BuildContext build_ctx_data;
173
174 build_ctx_data.dim = build_ctx_data.space_dim = dim;
175 CeedQFunctionContextCreate(ceed, &build_ctx);
176 CeedQFunctionContextSetData(build_ctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(build_ctx_data), &build_ctx_data);
177
178 // Create the QFunction that builds the diffusion operator (i.e. computes its quadrature data) and set its context data.
179 CeedQFunction qf_build;
180
181 if (gallery) {
182 // This creates the QFunction via the gallery.
183 char name[16] = "";
184 snprintf(name, sizeof name, "Poisson%" CeedInt_FMT "DBuild", dim);
185 CeedQFunctionCreateInteriorByName(ceed, name, &qf_build);
186 } else {
187 // This creates the QFunction directly.
188 CeedQFunctionCreateInterior(ceed, 1, build_diff, build_diff_loc, &qf_build);
189 CeedQFunctionAddInput(qf_build, "dx", num_comp_x * dim, CEED_EVAL_GRAD);
190 CeedQFunctionAddInput(qf_build, "weights", 1, CEED_EVAL_WEIGHT);
191 CeedQFunctionAddOutput(qf_build, "qdata", dim * (dim + 1) / 2, CEED_EVAL_NONE);
192 CeedQFunctionSetContext(qf_build, build_ctx);
193 }
194
195 // Create the operator that builds the quadrature data for the diffusion operator.
196 CeedOperator op_build;
197
198 CeedOperatorCreate(ceed, qf_build, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_build);
199 CeedOperatorSetField(op_build, "dx", mesh_restriction, mesh_basis, CEED_VECTOR_ACTIVE);
200 CeedOperatorSetField(op_build, "weights", CEED_ELEMRESTRICTION_NONE, mesh_basis, CEED_VECTOR_NONE);
201 CeedOperatorSetField(op_build, "qdata", q_data_restriction, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE);
202
203 // Compute the quadrature data for the diffusion operator.
204 CeedVector q_data;
205 CeedInt elem_qpts = CeedIntPow(num_qpts, dim);
206 CeedInt num_elem = 1;
207
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
215 if (gallery) {
216 // This creates the QFunction via the gallery.
217 char name[25] = "";
218 snprintf(name, sizeof name, "Poisson%" CeedInt_FMT "DApply", dim);
219 CeedQFunctionCreateInteriorByName(ceed, name, &qf_apply);
220 } else {
221 // This creates the QFunction directly.
222 CeedQFunctionCreateInterior(ceed, 1, apply_diff, apply_diff_loc, &qf_apply);
223 CeedQFunctionAddInput(qf_apply, "du", dim, CEED_EVAL_GRAD);
224 CeedQFunctionAddInput(qf_apply, "qdata", dim * (dim + 1) / 2, CEED_EVAL_NONE);
225 CeedQFunctionAddOutput(qf_apply, "dv", dim, CEED_EVAL_GRAD);
226 CeedQFunctionSetContext(qf_apply, build_ctx);
227 }
228
229 // Create the diffusion operator.
230 CeedOperator op_apply;
231
232 CeedOperatorCreate(ceed, qf_apply, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_apply);
233 CeedOperatorSetField(op_apply, "du", sol_restriction, sol_basis, CEED_VECTOR_ACTIVE);
234 CeedOperatorSetField(op_apply, "qdata", q_data_restriction, CEED_BASIS_NONE, q_data);
235 CeedOperatorSetField(op_apply, "dv", sol_restriction, sol_basis, CEED_VECTOR_ACTIVE);
236
237 // Create auxiliary solution-size vectors.
238 CeedVector u, v;
239
240 CeedVectorCreate(ceed, sol_size, &u);
241 CeedVectorCreate(ceed, sol_size, &v);
242
243 // Initialize 'u' with sum of coordinates, x+y+z.
244 {
245 CeedScalar *u_array;
246 const CeedScalar *x_array;
247
248 CeedVectorGetArrayWrite(u, CEED_MEM_HOST, &u_array);
249 CeedVectorGetArrayRead(mesh_coords, CEED_MEM_HOST, &x_array);
250 for (CeedInt i = 0; i < sol_size; i++) {
251 u_array[i] = 0;
252 for (CeedInt d = 0; d < dim; d++) u_array[i] += x_array[i + d * sol_size];
253 }
254 CeedVectorRestoreArray(u, &u_array);
255 CeedVectorRestoreArrayRead(mesh_coords, &x_array);
256 }
257
258 // Compute the mesh surface area using the diff operator: surface_area = 1^T \cdot abs( K \cdot x).
259 CeedOperatorApply(op_apply, u, v, CEED_REQUEST_IMMEDIATE);
260
261 // Benchmark runs
262 if (!test && benchmark) {
263 // LCOV_EXCL_START
264 printf(" Executing %d benchmarking runs...\n", benchmark);
265 // LCOV_EXCL_STOP
266 }
267 for (CeedInt i = 0; i < benchmark; i++) {
268 // LCOV_EXCL_START
269 CeedOperatorApply(op_apply, u, v, CEED_REQUEST_IMMEDIATE);
270 // LCOV_EXCL_STOP
271 }
272
273 // Compute and print the sum of the entries of 'v' giving the mesh surface area.
274 CeedScalar surface_area = 0.;
275 {
276 const CeedScalar *v_array;
277
278 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array);
279 for (CeedInt i = 0; i < sol_size; i++) surface_area += fabs(v_array[i]);
280 CeedVectorRestoreArrayRead(v, &v_array);
281 }
282 if (!test) {
283 // LCOV_EXCL_START
284 printf(" done.\n");
285 printf("Exact mesh surface area : % .14g\n", exact_surface_area);
286 printf("Computed mesh surface area : % .14g\n", surface_area);
287 printf("Surface area error : % .14g\n", surface_area - exact_surface_area);
288 // LCOV_EXCL_STOP
289 } else {
290 CeedScalar tol = (dim == 1 ? 10000. * CEED_EPSILON : dim == 2 ? 1E-1 : 1E-1);
291
292 if (fabs(surface_area - exact_surface_area) > tol) printf("Surface area error : % .14g\n", surface_area - exact_surface_area);
293 }
294
295 // Free dynamically allocated memory.
296 CeedVectorDestroy(&u);
297 CeedVectorDestroy(&v);
298 CeedVectorDestroy(&q_data);
299 CeedVectorDestroy(&mesh_coords);
300 CeedOperatorDestroy(&op_apply);
301 CeedQFunctionDestroy(&qf_apply);
302 CeedQFunctionContextDestroy(&build_ctx);
303 CeedOperatorDestroy(&op_build);
304 CeedQFunctionDestroy(&qf_build);
305 CeedElemRestrictionDestroy(&sol_restriction);
306 CeedElemRestrictionDestroy(&mesh_restriction);
307 CeedElemRestrictionDestroy(&q_data_restriction);
308 CeedBasisDestroy(&sol_basis);
309 CeedBasisDestroy(&mesh_basis);
310 CeedDestroy(&ceed);
311 return 0;
312 }
313
GetCartesianMeshSize(CeedInt dim,CeedInt degree,CeedInt prob_size,CeedInt num_xyz[3])314 int GetCartesianMeshSize(CeedInt dim, CeedInt degree, CeedInt prob_size, CeedInt num_xyz[3]) {
315 // Use the approximate formula:
316 // prob_size ~ num_elem * degree^dim
317 CeedInt num_elem = prob_size / CeedIntPow(degree, dim);
318 CeedInt s = 0; // find s: num_elem/2 < 2^s <= num_elem
319
320 while (num_elem > 1) {
321 num_elem /= 2;
322 s++;
323 }
324 CeedInt r = s % dim;
325
326 for (CeedInt d = 0; d < dim; d++) {
327 CeedInt sd = s / dim;
328
329 if (r > 0) {
330 sd++;
331 r--;
332 }
333 num_xyz[d] = 1 << sd;
334 }
335 return 0;
336 }
337
BuildCartesianRestriction(Ceed ceed,CeedInt dim,CeedInt num_xyz[3],CeedInt degree,CeedInt num_comp,CeedInt * size,CeedInt num_qpts,CeedElemRestriction * restriction,CeedElemRestriction * q_data_restriction)338 int BuildCartesianRestriction(Ceed ceed, CeedInt dim, CeedInt num_xyz[3], CeedInt degree, CeedInt num_comp, CeedInt *size, CeedInt num_qpts,
339 CeedElemRestriction *restriction, CeedElemRestriction *q_data_restriction) {
340 CeedInt p = degree + 1;
341 CeedInt num_nodes = CeedIntPow(p, dim); // number of scalar nodes per element
342 CeedInt elem_qpts = CeedIntPow(num_qpts, dim); // number of qpts per element
343 CeedInt nd[3], num_elem = 1, scalar_size = 1;
344
345 for (CeedInt d = 0; d < dim; d++) {
346 num_elem *= num_xyz[d];
347 nd[d] = num_xyz[d] * (p - 1) + 1;
348 scalar_size *= nd[d];
349 }
350 *size = scalar_size * num_comp;
351 // elem: 0 1 n-1
352 // |---*-...-*---|---*-...-*---|- ... -|--...--|
353 // num_nodes: 0 1 p-1 p p+1 2*p n*p
354 CeedInt *el_nodes = malloc(sizeof(CeedInt) * num_elem * num_nodes);
355
356 for (CeedInt e = 0; e < num_elem; e++) {
357 CeedInt e_xyz[3] = {1, 1, 1}, re = e;
358
359 for (CeedInt d = 0; d < dim; d++) {
360 e_xyz[d] = re % num_xyz[d];
361 re /= num_xyz[d];
362 }
363 CeedInt *local_elem_nodes = el_nodes + e * num_nodes;
364
365 for (CeedInt l_nodes = 0; l_nodes < num_nodes; l_nodes++) {
366 CeedInt g_nodes = 0, g_nodes_stride = 1, r_nodes = l_nodes;
367
368 for (CeedInt d = 0; d < dim; d++) {
369 g_nodes += (e_xyz[d] * (p - 1) + r_nodes % p) * g_nodes_stride;
370 g_nodes_stride *= nd[d];
371 r_nodes /= p;
372 }
373 local_elem_nodes[l_nodes] = g_nodes;
374 }
375 }
376 if (restriction) {
377 CeedElemRestrictionCreate(ceed, num_elem, num_nodes, num_comp, scalar_size, num_comp * scalar_size, CEED_MEM_HOST, CEED_COPY_VALUES, el_nodes,
378 restriction);
379 }
380 free(el_nodes);
381
382 if (q_data_restriction) {
383 CeedElemRestrictionCreateStrided(ceed, num_elem, elem_qpts, num_comp, num_comp * elem_qpts * num_elem, CEED_STRIDES_BACKEND, q_data_restriction);
384 }
385
386 return 0;
387 }
388
SetCartesianMeshCoords(CeedInt dim,CeedInt num_xyz[3],CeedInt mesh_degree,CeedVector mesh_coords)389 int SetCartesianMeshCoords(CeedInt dim, CeedInt num_xyz[3], CeedInt mesh_degree, CeedVector mesh_coords) {
390 CeedInt p = mesh_degree + 1;
391 CeedInt nd[3], scalar_size = 1;
392
393 for (CeedInt d = 0; d < dim; d++) {
394 nd[d] = num_xyz[d] * (p - 1) + 1;
395 scalar_size *= nd[d];
396 }
397 CeedScalar *coords;
398
399 CeedVectorGetArrayWrite(mesh_coords, CEED_MEM_HOST, &coords);
400 CeedScalar *nodes = malloc(sizeof(CeedScalar) * p);
401
402 // The H1 basis uses Lobatto quadrature points as nodes.
403 CeedLobattoQuadrature(p, nodes, NULL); // nodes are in [-1,1]
404 for (CeedInt i = 0; i < p; i++) nodes[i] = 0.5 + 0.5 * nodes[i];
405 for (CeedInt gs_nodes = 0; gs_nodes < scalar_size; gs_nodes++) {
406 CeedInt r_nodes = gs_nodes;
407
408 for (CeedInt d = 0; d < dim; d++) {
409 CeedInt d1d = r_nodes % nd[d];
410
411 coords[gs_nodes + scalar_size * d] = ((d1d / (p - 1)) + nodes[d1d % (p - 1)]) / num_xyz[d];
412 r_nodes /= nd[d];
413 }
414 }
415 free(nodes);
416 CeedVectorRestoreArray(mesh_coords, &coords);
417 return 0;
418 }
419
420 #ifndef M_PI
421 #define M_PI 3.14159265358979323846
422 #endif
423
TransformMeshCoords(CeedInt dim,CeedInt mesh_size,CeedVector mesh_coords)424 CeedScalar TransformMeshCoords(CeedInt dim, CeedInt mesh_size, CeedVector mesh_coords) {
425 CeedScalar exact_surface_area = (dim == 1 ? 2 : dim == 2 ? 4 : 6);
426 CeedScalar *coords;
427
428 CeedVectorGetArray(mesh_coords, CEED_MEM_HOST, &coords);
429 for (CeedInt i = 0; i < mesh_size; i++) {
430 // map [0,1] to [0,1] varying the mesh density
431 coords[i] = 0.5 + 1. / sqrt(3.) * sin((2. / 3.) * M_PI * (coords[i] - 0.5));
432 }
433 CeedVectorRestoreArray(mesh_coords, &coords);
434 return exact_surface_area;
435 }
436