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