1 // Copyright (c) 2017-2024, 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 operator. Arbitrary mesh and 12 // 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::{prelude::*, Ceed}; 23 mod opt; 24 mod transform; 25 26 // ---------------------------------------------------------------------------- 27 // Example 1 28 // ---------------------------------------------------------------------------- 29 fn main() -> libceed::Result<()> { 30 let options = opt::Opt::parse(); 31 example_1(options) 32 } 33 34 fn example_1(options: opt::Opt) -> libceed::Result<()> { 35 // Process command line arguments 36 let opt::Opt { 37 ceed_spec, 38 dim, 39 mesh_degree, 40 solution_degree, 41 num_qpts, 42 problem_size_requested, 43 test, 44 quiet, 45 gallery, 46 } = options; 47 assert!(dim >= 1 && dim <= 3); 48 assert!(mesh_degree >= 1); 49 assert!(solution_degree >= 1); 50 assert!(num_qpts >= 1); 51 let ncomp_x = dim; 52 let problem_size: i64; 53 if problem_size_requested < 0 { 54 problem_size = if test { 8 * 16 } else { 256 * 1024 }; 55 } else { 56 problem_size = problem_size_requested; 57 } 58 59 // Summary output 60 if !quiet { 61 println!("Selected options: [command line option] : <current value>"); 62 println!(" Ceed specification [-c] : {}", ceed_spec); 63 println!(" Mesh dimension [-d] : {}", dim); 64 println!(" Mesh degree [-m] : {}", mesh_degree); 65 println!(" Solution degree [-p] : {}", solution_degree); 66 println!(" Num. 1D quadr. pts [-q] : {}", num_qpts); 67 println!(" Approx. # unknowns [-s] : {}", problem_size); 68 println!( 69 " QFunction source [-g] : {}", 70 if gallery { "gallery" } else { "user closure" } 71 ); 72 } 73 74 // Initalize ceed context 75 let ceed = Ceed::init(&ceed_spec); 76 77 // Mesh and solution bases 78 let basis_mesh = 79 ceed.basis_tensor_H1_Lagrange(dim, ncomp_x, mesh_degree + 1, num_qpts, QuadMode::Gauss)?; 80 let basis_solution = 81 ceed.basis_tensor_H1_Lagrange(dim, 1, solution_degree + 1, num_qpts, QuadMode::Gauss)?; 82 83 // Determine mesh size from approximate problem size 84 let num_xyz = mesh::cartesian_mesh_size(dim, solution_degree, problem_size); 85 if !quiet { 86 print!("\nMesh size : nx = {}", num_xyz[0]); 87 if dim > 1 { 88 print!(", ny = {}", num_xyz[1]); 89 } 90 if dim > 2 { 91 print!(", nz = {}", num_xyz[2]); 92 } 93 print!("\n"); 94 } 95 96 // Build ElemRestriction objects describing the mesh and solution discrete 97 // representations 98 let (rstr_mesh, _) = 99 mesh::build_cartesian_restriction(&ceed, dim, num_xyz, mesh_degree, ncomp_x, num_qpts)?; 100 let (rstr_solution, rstr_qdata) = 101 mesh::build_cartesian_restriction(&ceed, dim, num_xyz, solution_degree, 1, num_qpts)?; 102 let mesh_size = rstr_mesh.lvector_size(); 103 let solution_size = rstr_solution.lvector_size(); 104 if !quiet { 105 println!("Number of mesh nodes : {}", mesh_size / dim); 106 println!("Number of solution nodes : {}", solution_size); 107 } 108 109 // Create a Vector with the mesh coordinates 110 let mut mesh_coords = mesh::cartesian_mesh_coords(&ceed, dim, num_xyz, mesh_degree, mesh_size)?; 111 112 // Apply a transformation to the mesh coordinates 113 let exact_volume = transform::transform_mesh_coordinates(dim, mesh_size, &mut mesh_coords)?; 114 115 // QFunction that builds the quadrature data for the mass operator 116 // -- QFunction from user closure 117 let build_mass = move |[jacobian, weights, ..]: QFunctionInputs, 118 [qdata, ..]: QFunctionOutputs| { 119 // Build quadrature data 120 match dim { 121 1 => qdata 122 .iter_mut() 123 .zip(jacobian.iter().zip(weights.iter())) 124 .for_each(|(qdata, (j, weight))| *qdata = j * weight), 125 2 => { 126 let q = qdata.len(); 127 qdata.iter_mut().zip(weights.iter()).enumerate().for_each( 128 |(i, (qdata, weight))| { 129 *qdata = (jacobian[i + q * 0] * jacobian[i + q * 3] 130 - jacobian[i + q * 1] * jacobian[i + q * 2]) 131 * weight 132 }, 133 ); 134 } 135 3 => { 136 let q = qdata.len(); 137 qdata.iter_mut().zip(weights.iter()).enumerate().for_each( 138 |(i, (qdata, weight))| { 139 *qdata = (jacobian[i + q * 0] 140 * (jacobian[i + q * 4] * jacobian[i + q * 8] 141 - jacobian[i + q * 5] * jacobian[i + q * 7]) 142 - jacobian[i + q * 1] 143 * (jacobian[i + q * 3] * jacobian[i + q * 8] 144 - jacobian[i + q * 5] * jacobian[i + q * 6]) 145 + jacobian[i + q * 2] 146 * (jacobian[i + q * 3] * jacobian[i + q * 7] 147 - jacobian[i + q * 4] * jacobian[i + q * 6])) 148 * weight 149 }, 150 ); 151 } 152 _ => unreachable!(), 153 }; 154 155 // Return clean error code 156 0 157 }; 158 let qf_build_closure = ceed 159 .q_function_interior(1, Box::new(build_mass))? 160 .input("dx", ncomp_x * dim, EvalMode::Grad)? 161 .input("weights", 1, EvalMode::Weight)? 162 .output("qdata", 1, EvalMode::None)?; 163 // -- QFunction from gallery 164 let qf_build_named = { 165 let name = format!("Mass{}DBuild", dim); 166 ceed.q_function_interior_by_name(&name)? 167 }; 168 // -- QFunction for use with Operator 169 let qf_build = if gallery { 170 QFunctionOpt::SomeQFunctionByName(&qf_build_named) 171 } else { 172 QFunctionOpt::SomeQFunction(&qf_build_closure) 173 }; 174 175 // Operator that build the quadrature data for the mass operator 176 let op_build = ceed 177 .operator(qf_build, QFunctionOpt::None, QFunctionOpt::None)? 178 .name("build qdata")? 179 .field("dx", &rstr_mesh, &basis_mesh, VectorOpt::Active)? 180 .field( 181 "weights", 182 ElemRestrictionOpt::None, 183 &basis_mesh, 184 VectorOpt::None, 185 )? 186 .field("qdata", &rstr_qdata, BasisOpt::None, VectorOpt::Active)? 187 .check()?; 188 189 // Compute the quadrature data for the mass operator 190 let elem_qpts = num_qpts.pow(dim as u32); 191 let num_elem: usize = num_xyz.iter().take(dim).product(); 192 let mut qdata = ceed.vector(num_elem * elem_qpts)?; 193 op_build.apply(&mesh_coords, &mut qdata)?; 194 195 // QFunction that applies the mass operator 196 // -- QFunction from user closure 197 let apply_mass = |[u, qdata, ..]: QFunctionInputs, [v, ..]: QFunctionOutputs| { 198 // Apply mass operator 199 v.iter_mut() 200 .zip(u.iter().zip(qdata.iter())) 201 .for_each(|(v, (u, w))| *v = u * w); 202 // Return clean error code 203 0 204 }; 205 let qf_mass_closure = ceed 206 .q_function_interior(1, Box::new(apply_mass))? 207 .input("u", 1, EvalMode::Interp)? 208 .input("qdata", 1, EvalMode::None)? 209 .output("v", 1, EvalMode::Interp)?; 210 // -- QFunction from gallery 211 let qf_mass_named = ceed.q_function_interior_by_name("MassApply")?; 212 // -- QFunction for use with Operator 213 let qf_mass = if gallery { 214 QFunctionOpt::SomeQFunctionByName(&qf_mass_named) 215 } else { 216 QFunctionOpt::SomeQFunction(&qf_mass_closure) 217 }; 218 219 // Mass Operator 220 let op_mass = ceed 221 .operator(qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 222 .name("mass")? 223 .field("u", &rstr_solution, &basis_solution, VectorOpt::Active)? 224 .field("qdata", &rstr_qdata, BasisOpt::None, &qdata)? 225 .field("v", &rstr_solution, &basis_solution, VectorOpt::Active)? 226 .check()?; 227 228 // Solution vectors 229 let u = ceed.vector_from_slice(&vec![1.0; solution_size])?; 230 let mut v = ceed.vector(solution_size)?; 231 232 // Apply the mass operator 233 op_mass.apply(&u, &mut v)?; 234 235 // Compute the mesh volume 236 let volume: Scalar = v.view()?.iter().sum(); 237 238 // Output results 239 if !quiet { 240 println!("Exact mesh volume : {:.12}", exact_volume); 241 println!("Computed mesh volume : {:.12}", volume); 242 println!( 243 "Volume error : {:.12e}", 244 volume - exact_volume 245 ); 246 } 247 let tolerance = match dim { 248 1 => 500.0 * libceed::EPSILON, 249 _ => 1E-5, 250 }; 251 let error = (volume - exact_volume).abs(); 252 if error > tolerance { 253 println!("Volume error too large: {:.12e}", error); 254 return Err(libceed::Error { 255 message: format!( 256 "Volume error too large - expected: {:.12e}, actual: {:.12e}", 257 tolerance, error 258 ), 259 }); 260 } 261 Ok(()) 262 } 263 264 // ---------------------------------------------------------------------------- 265 // Tests 266 // ---------------------------------------------------------------------------- 267 #[cfg(test)] 268 mod tests { 269 use super::*; 270 271 #[test] 272 fn example_1_1d() { 273 let options = opt::Opt { 274 ceed_spec: "/cpu/self/ref/serial".to_string(), 275 dim: 1, 276 mesh_degree: 4, 277 solution_degree: 4, 278 num_qpts: 6, 279 problem_size_requested: -1, 280 test: true, 281 quiet: true, 282 gallery: false, 283 }; 284 assert!(example_1(options).is_ok()); 285 } 286 287 #[test] 288 fn example_1_2d() { 289 let options = opt::Opt { 290 ceed_spec: "/cpu/self/ref/serial".to_string(), 291 dim: 2, 292 mesh_degree: 4, 293 solution_degree: 4, 294 num_qpts: 6, 295 problem_size_requested: -1, 296 test: true, 297 quiet: true, 298 gallery: false, 299 }; 300 assert!(example_1(options).is_ok()); 301 } 302 303 #[test] 304 fn example_1_3d() { 305 let options = opt::Opt { 306 ceed_spec: "/cpu/self/ref/serial".to_string(), 307 dim: 3, 308 mesh_degree: 4, 309 solution_degree: 4, 310 num_qpts: 6, 311 problem_size_requested: -1, 312 test: true, 313 quiet: false, 314 gallery: false, 315 }; 316 assert!(example_1(options).is_ok()); 317 } 318 319 #[test] 320 fn example_1_1d_gallery() { 321 let options = opt::Opt { 322 ceed_spec: "/cpu/self/ref/serial".to_string(), 323 dim: 1, 324 mesh_degree: 4, 325 solution_degree: 4, 326 num_qpts: 6, 327 problem_size_requested: -1, 328 test: true, 329 quiet: true, 330 gallery: true, 331 }; 332 assert!(example_1(options).is_ok()); 333 } 334 335 #[test] 336 fn example_1_2d_gallery() { 337 let options = opt::Opt { 338 ceed_spec: "/cpu/self/ref/serial".to_string(), 339 dim: 2, 340 mesh_degree: 4, 341 solution_degree: 4, 342 num_qpts: 6, 343 problem_size_requested: -1, 344 test: true, 345 quiet: true, 346 gallery: true, 347 }; 348 assert!(example_1(options).is_ok()); 349 } 350 351 #[test] 352 fn example_1_3d_gallery() { 353 let options = opt::Opt { 354 ceed_spec: "/cpu/self/ref/serial".to_string(), 355 dim: 3, 356 mesh_degree: 4, 357 solution_degree: 4, 358 num_qpts: 6, 359 problem_size_requested: -1, 360 test: true, 361 quiet: true, 362 gallery: true, 363 }; 364 assert!(example_1(options).is_ok()); 365 } 366 } 367 368 // ---------------------------------------------------------------------------- 369