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