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