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