1 // Copyright (c) 2017-2025, 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 // ---------------------------------------------------------------------------- 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)] 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] 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] 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] 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