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
11 // 3D body using matrix-free application of a mass + diff operator. Arbitrary
12 // mesh and solution orders in 1D, 2D and 3D are supported from the same code.
13 // This calculation is executed in triplicate with a 3 component vector system.
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 use clap::Parser;
23 use libceed::{
24 BasisOpt, Ceed, ElemRestrictionOpt, QFunctionInputs, QFunctionOpt, QFunctionOutputs, VectorOpt,
25 };
26 mod opt;
27 mod transform;
28
29 // ----------------------------------------------------------------------------
30 // Example 1
31 // ----------------------------------------------------------------------------
main() -> libceed::Result<()>32 fn main() -> libceed::Result<()> {
33 let options = opt::Opt::parse();
34 example_3_vector(options)
35 }
36
37 #[allow(clippy::erasing_op)]
38 #[allow(clippy::identity_op)]
example_3_vector(options: opt::Opt) -> libceed::Result<()>39 fn example_3_vector(options: opt::Opt) -> libceed::Result<()> {
40 // Process command line arguments
41 let opt::Opt {
42 ceed_spec,
43 dim,
44 mesh_degree,
45 solution_degree,
46 num_qpts,
47 problem_size_requested,
48 test,
49 quiet,
50 } = options;
51 assert!((1..=3).contains(&dim));
52 assert!(mesh_degree >= 1);
53 assert!(solution_degree >= 1);
54 assert!(num_qpts >= 1);
55 let ncomp_x = dim;
56 let problem_size: i64 = if problem_size_requested < 0 {
57 if test {
58 8 * 16
59 } else {
60 256 * 1024
61 }
62 } else {
63 problem_size_requested
64 };
65 let ncomp_u = 3;
66
67 // Summary output
68 if !quiet {
69 println!("Selected options: [command line option] : <current value>");
70 println!(" Ceed specification [-c] : {}", ceed_spec);
71 println!(" Mesh dimension [-d] : {}", dim);
72 println!(" Mesh degree [-m] : {}", mesh_degree);
73 println!(" Solution degree [-p] : {}", solution_degree);
74 println!(" Num. 1D quadr. pts [-q] : {}", num_qpts);
75 println!(" Approx. # unknowns [-s] : {}", problem_size);
76 println!(" QFunction source : user closure");
77 }
78
79 // Initalize ceed context
80 let ceed = Ceed::init(&ceed_spec);
81
82 // Mesh and solution bases
83 let basis_mesh = ceed.basis_tensor_H1_Lagrange(
84 dim,
85 ncomp_x,
86 mesh_degree + 1,
87 num_qpts,
88 libceed::QuadMode::Gauss,
89 )?;
90 let basis_solution = ceed.basis_tensor_H1_Lagrange(
91 dim,
92 ncomp_u,
93 solution_degree + 1,
94 num_qpts,
95 libceed::QuadMode::Gauss,
96 )?;
97
98 // Determine mesh size from approximate problem size
99 let num_xyz = mesh::cartesian_mesh_size(dim, solution_degree, problem_size);
100 if !quiet {
101 print!("\nMesh size : nx = {}", num_xyz[0]);
102 if dim > 1 {
103 print!(", ny = {}", num_xyz[1]);
104 }
105 if dim > 2 {
106 print!(", nz = {}", num_xyz[2]);
107 }
108 println!();
109 }
110
111 // Build ElemRestriction objects describing the mesh and solution discrete
112 // representations
113 let (rstr_mesh, _) =
114 mesh::build_cartesian_restriction(&ceed, dim, num_xyz, mesh_degree, ncomp_x, num_qpts)?;
115 let (_, rstr_qdata) = mesh::build_cartesian_restriction(
116 &ceed,
117 dim,
118 num_xyz,
119 solution_degree,
120 1 + dim * (dim + 1) / 2,
121 num_qpts,
122 )?;
123 let (rstr_solution, _) =
124 mesh::build_cartesian_restriction(&ceed, dim, num_xyz, solution_degree, ncomp_u, num_qpts)?;
125 let mesh_size = rstr_mesh.lvector_size();
126 let solution_size = rstr_solution.lvector_size();
127 if !quiet {
128 println!("Number of mesh nodes : {}", mesh_size / dim);
129 println!("Number of solution nodes : {}", solution_size);
130 }
131
132 // Create a Vector with the mesh coordinates
133 let mut mesh_coords = mesh::cartesian_mesh_coords(&ceed, dim, num_xyz, mesh_degree, mesh_size)?;
134
135 // Apply a transformation to the mesh coordinates
136 let exact_volume = transform::transform_mesh_coordinates(dim, mesh_size, &mut mesh_coords)?;
137
138 // QFunction that builds the quadrature data for the mass + diff operator
139 // -- QFunction from user closure
140 let build_mass_diff = move |[jacobian, weights, ..]: QFunctionInputs,
141 [qdata, ..]: QFunctionOutputs| {
142 // Build quadrature data
143 match dim {
144 1 => {
145 let q = qdata.len() / 2;
146 for i in 0..q {
147 // Mass
148 qdata[i + q * 0] = weights[i] * jacobian[i];
149 // Diff
150 qdata[i + q * 1] = weights[i] / jacobian[i];
151 }
152 }
153 2 => {
154 let q = qdata.len() / 4;
155 for i in 0..q {
156 let j11 = jacobian[i + q * 0];
157 let j21 = jacobian[i + q * 1];
158 let j12 = jacobian[i + q * 2];
159 let j22 = jacobian[i + q * 3];
160 // Mass
161 qdata[i + q * 0] = weights[i] * (j11 * j22 - j21 * j12);
162 // Diff
163 let qw = weights[i] / (j11 * j22 - j21 * j12);
164 qdata[i + q * 1] = qw * (j12 * j12 + j22 * j22);
165 qdata[i + q * 2] = qw * (j11 * j11 + j21 * j21);
166 qdata[i + q * 3] = -qw * (j11 * j12 + j21 * j22);
167 }
168 }
169 3 => {
170 let q = qdata.len() / 7;
171 for i in 0..q {
172 let mut a = [0.0; 9];
173 for j in 0..3 {
174 for k in 0..3 {
175 a[k * 3 + j] = jacobian[i + q * ((j + 1) % 3 + 3 * ((k + 1) % 3))]
176 * jacobian[i + q * ((j + 2) % 3 + 3 * ((k + 2) % 3))]
177 - jacobian[i + q * ((j + 1) % 3 + 3 * ((k + 2) % 3))]
178 * jacobian[i + q * ((j + 2) % 3 + 3 * ((k + 1) % 3))];
179 }
180 }
181 // Mass
182 qdata[i + q * 0] = weights[i]
183 * (jacobian[i + q * 0] * a[0 * 3 + 0]
184 + jacobian[i + q * 1] * a[0 * 3 + 1]
185 + jacobian[i + q * 2] * a[0 * 3 + 2]);
186 let qw = weights[i]
187 / (jacobian[i + q * 0] * a[0 * 3 + 0]
188 + jacobian[i + q * 1] * a[0 * 3 + 1]
189 + jacobian[i + q * 2] * a[0 * 3 + 2]);
190 // Diff
191 qdata[i + q * 1] = qw
192 * (a[0 * 3 + 0] * a[0 * 3 + 0]
193 + a[0 * 3 + 1] * a[0 * 3 + 1]
194 + a[0 * 3 + 2] * a[0 * 3 + 2]);
195 qdata[i + q * 2] = qw
196 * (a[1 * 3 + 0] * a[1 * 3 + 0]
197 + a[1 * 3 + 1] * a[1 * 3 + 1]
198 + a[1 * 3 + 2] * a[1 * 3 + 2]);
199 qdata[i + q * 3] = qw
200 * (a[2 * 3 + 0] * a[2 * 3 + 0]
201 + a[2 * 3 + 1] * a[2 * 3 + 1]
202 + a[2 * 3 + 2] * a[2 * 3 + 2]);
203 qdata[i + q * 4] = qw
204 * (a[1 * 3 + 0] * a[2 * 3 + 0]
205 + a[1 * 3 + 1] * a[2 * 3 + 1]
206 + a[1 * 3 + 2] * a[2 * 3 + 2]);
207 qdata[i + q * 5] = qw
208 * (a[0 * 3 + 0] * a[2 * 3 + 0]
209 + a[0 * 3 + 1] * a[2 * 3 + 1]
210 + a[0 * 3 + 2] * a[2 * 3 + 2]);
211 qdata[i + q * 6] = qw
212 * (a[0 * 3 + 0] * a[1 * 3 + 0]
213 + a[0 * 3 + 1] * a[1 * 3 + 1]
214 + a[0 * 3 + 2] * a[1 * 3 + 2]);
215 }
216 }
217 _ => unreachable!(),
218 };
219
220 // Return clean error code
221 0
222 };
223 let qf_build_closure = ceed
224 .q_function_interior(1, Box::new(build_mass_diff))?
225 .input("dx", ncomp_x * dim, libceed::EvalMode::Grad)?
226 .input("weights", 1, libceed::EvalMode::Weight)?
227 .output("qdata", 1 + dim * (dim + 1) / 2, libceed::EvalMode::None)?;
228 // -- QFunction for use with Operator
229 let qf_build = QFunctionOpt::SomeQFunction(&qf_build_closure);
230
231 // Operator that build the quadrature data for the mass + diff operator
232 let op_build = ceed
233 .operator(qf_build, QFunctionOpt::None, QFunctionOpt::None)?
234 .name("build qdata")?
235 .field("dx", &rstr_mesh, &basis_mesh, VectorOpt::Active)?
236 .field(
237 "weights",
238 ElemRestrictionOpt::None,
239 &basis_mesh,
240 VectorOpt::None,
241 )?
242 .field("qdata", &rstr_qdata, BasisOpt::None, VectorOpt::Active)?
243 .check()?;
244
245 // Compute the quadrature data for the mass + diff operator
246 let elem_qpts = num_qpts.pow(dim as u32);
247 let num_elem: usize = num_xyz.iter().take(dim).product();
248 let mut qdata = ceed.vector(num_elem * elem_qpts * (1 + dim * (dim + 1) / 2))?;
249 op_build.apply(&mesh_coords, &mut qdata)?;
250
251 // QFunction that applies the mass + diff operator
252 // -- QFunction from user closure
253 let apply_mass_diff = move |[u, ug, qdata, ..]: QFunctionInputs,
254 [v, vg, ..]: QFunctionOutputs| {
255 // Apply diffusion operator
256 match dim {
257 1 => {
258 let q = qdata.len() / 2;
259 for i in 0..q {
260 for c in 0..ncomp_u {
261 // Mass
262 v[i + c * q] = u[i + c * q] * qdata[i + 0 * q];
263 // Diff
264 vg[i + c * q] = ug[i + c * q] * qdata[i + 1 * q];
265 }
266 }
267 }
268 2 => {
269 let q = qdata.len() / 4;
270 for i in 0..q {
271 let dxdxdxdx_t = [
272 [qdata[i + 1 * q], qdata[i + 3 * q]],
273 [qdata[i + 3 * q], qdata[i + 2 * q]],
274 ];
275 for c in 0..ncomp_u {
276 // Mass
277 v[i + c * q] = u[i + c * q] * qdata[i + 0 * q];
278 // Diff
279 let du = [ug[i + (c + 0 * ncomp_u) * q], ug[i + (c + 1 * ncomp_u) * q]];
280 for j in 0..2 {
281 vg[i + (j + j * ncomp_u) * q] =
282 du[0] * dxdxdxdx_t[0][j] + du[1] * dxdxdxdx_t[1][j];
283 }
284 }
285 }
286 }
287 3 => {
288 let q = qdata.len() / 7;
289 for i in 0..q {
290 let dxdxdxdx_t = [
291 [qdata[i + 1 * q], qdata[i + 6 * q], qdata[i + 5 * q]],
292 [qdata[i + 6 * q], qdata[i + 2 * q], qdata[i + 4 * q]],
293 [qdata[i + 5 * q], qdata[i + 4 * q], qdata[i + 3 * q]],
294 ];
295 for c in 0..ncomp_u {
296 // Mass
297 v[i + c * q] = u[i + c * q] * qdata[i + 0 * q];
298 // Diff
299 let du = [
300 ug[i + (c + 0 * ncomp_u) * q],
301 ug[i + (c + 1 * ncomp_u) * q],
302 ug[i + (c + 2 * ncomp_u) * q],
303 ];
304 for j in 0..3 {
305 vg[i + (c + j * ncomp_u) * q] = du[0] * dxdxdxdx_t[0][j]
306 + du[1] * dxdxdxdx_t[1][j]
307 + du[2] * dxdxdxdx_t[2][j];
308 }
309 }
310 }
311 }
312 _ => unreachable!(),
313 };
314
315 // Return clean error code
316 0
317 };
318 let qf_mass_diff_closure = ceed
319 .q_function_interior(1, Box::new(apply_mass_diff))?
320 .input("u", ncomp_u, libceed::EvalMode::Interp)?
321 .input("du", dim * ncomp_u, libceed::EvalMode::Grad)?
322 .input("qdata", 1 + dim * (dim + 1) / 2, libceed::EvalMode::None)?
323 .output("v", ncomp_u, libceed::EvalMode::Interp)?
324 .output("dv", dim * ncomp_u, libceed::EvalMode::Grad)?;
325 // -- QFunction for use with Operator
326 let qf_mass_diff = QFunctionOpt::SomeQFunction(&qf_mass_diff_closure);
327
328 // Mass + diff Operator
329 let op_mass_diff = ceed
330 .operator(qf_mass_diff, QFunctionOpt::None, QFunctionOpt::None)?
331 .name("mass diff")?
332 .field("u", &rstr_solution, &basis_solution, VectorOpt::Active)?
333 .field("du", &rstr_solution, &basis_solution, VectorOpt::Active)?
334 .field("qdata", &rstr_qdata, BasisOpt::None, &qdata)?
335 .field("v", &rstr_solution, &basis_solution, VectorOpt::Active)?
336 .field("dv", &rstr_solution, &basis_solution, VectorOpt::Active)?
337 .check()?;
338
339 // Solution vectors
340 let mut u = ceed.vector(solution_size)?;
341 let mut v = ceed.vector(solution_size)?;
342
343 // Initialize u with component index
344 u.set_value(0.0)?;
345 for c in 0..ncomp_u {
346 let q = solution_size / ncomp_u;
347 u.view_mut()?.iter_mut().skip(c * q).take(q).for_each(|u| {
348 *u = (c + 1) as libceed::Scalar;
349 });
350 }
351
352 // Apply the mass + diff operator
353 op_mass_diff.apply(&u, &mut v)?;
354
355 // Compute the mesh volume
356 let volume: libceed::Scalar = v.view()?.iter().sum::<libceed::Scalar>()
357 / ((ncomp_u * (ncomp_u + 1)) / 2) as libceed::Scalar;
358
359 // Output results
360 if !quiet {
361 println!("Exact mesh volume : {:.12}", exact_volume);
362 println!("Computed mesh volume : {:.12}", volume);
363 println!(
364 "Volume error : {:.12e}",
365 volume - exact_volume
366 );
367 }
368 let tolerance = match dim {
369 1 => 200.0 * libceed::EPSILON,
370 _ => 1E-5,
371 };
372 let error = (volume - exact_volume).abs();
373 if error > tolerance {
374 println!("Volume error too large: {:.12e}", error);
375 return Err(libceed::Error {
376 message: format!(
377 "Volume error too large - expected: {:.12e}, actual: {:.12e}",
378 tolerance, error
379 ),
380 });
381 }
382 Ok(())
383 }
384
385 // ----------------------------------------------------------------------------
386 // Tests
387 // ----------------------------------------------------------------------------
388 #[cfg(test)]
389 mod tests {
390 use super::*;
391
392 #[test]
example_3_vector_1d()393 fn example_3_vector_1d() {
394 let options = opt::Opt {
395 ceed_spec: "/cpu/self/ref/serial".to_string(),
396 dim: 1,
397 mesh_degree: 4,
398 solution_degree: 4,
399 num_qpts: 6,
400 problem_size_requested: -1,
401 test: true,
402 quiet: true,
403 };
404 assert!(example_3_vector(options).is_ok());
405 }
406
407 #[test]
example_3_vector_2d()408 fn example_3_vector_2d() {
409 let options = opt::Opt {
410 ceed_spec: "/cpu/self/ref/serial".to_string(),
411 dim: 2,
412 mesh_degree: 4,
413 solution_degree: 4,
414 num_qpts: 6,
415 problem_size_requested: -1,
416 test: true,
417 quiet: true,
418 };
419 assert!(example_3_vector(options).is_ok());
420 }
421
422 #[test]
example_3_vector_vector_3d()423 fn example_3_vector_vector_3d() {
424 let options = opt::Opt {
425 ceed_spec: "/cpu/self/ref/serial".to_string(),
426 dim: 3,
427 mesh_degree: 4,
428 solution_degree: 4,
429 num_qpts: 6,
430 problem_size_requested: -1,
431 test: true,
432 quiet: false,
433 };
434 assert!(example_3_vector(options).is_ok());
435 }
436 }
437
438 // ----------------------------------------------------------------------------
439