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