1 // Copyright (c) 2017-2021, Lawrence Livermore National Security, LLC. 2 // Produced at the Lawrence Livermore National Laboratory. LLNL-CODE-734707. 3 // All Rights reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 // libCEED Example 1 18 // 19 // This example illustrates a simple usage of libCEED to compute the volume of a 20 // 3D body using matrix-free application of a mass operator. Arbitrary mesh and 21 // solution orders in 1D, 2D and 3D are supported from the same code. 22 // 23 // The example has no dependencies, and is designed to be self-contained. For 24 // additional examples that use external discretization libraries (MFEM, PETSc, 25 // etc.) see the subdirectories in libceed/examples. 26 // 27 // All libCEED objects use a Ceed device object constructed based on a command 28 // line argument (-ceed). 29 30 use libceed::{prelude::*, Ceed}; 31 use mesh; 32 use structopt::StructOpt; 33 34 mod opt; 35 mod transform; 36 37 // ---------------------------------------------------------------------------- 38 // Example 1 39 // ---------------------------------------------------------------------------- 40 #[cfg(not(tarpaulin_include))] 41 fn main() -> libceed::Result<()> { 42 let options = opt::Opt::from_args(); 43 example_1(options) 44 } 45 46 fn example_1(options: opt::Opt) -> libceed::Result<()> { 47 // Process command line arguments 48 let opt::Opt { 49 ceed_spec, 50 dim, 51 mesh_degree, 52 solution_degree, 53 num_qpts, 54 problem_size_requested, 55 test, 56 quiet, 57 gallery, 58 } = options; 59 assert!(dim >= 1 && dim <= 3); 60 assert!(mesh_degree >= 1); 61 assert!(solution_degree >= 1); 62 assert!(num_qpts >= 1); 63 let ncomp_x = dim; 64 let problem_size: i64; 65 if problem_size_requested < 0 { 66 problem_size = if test { 8 * 16 } else { 256 * 1024 }; 67 } else { 68 problem_size = problem_size_requested; 69 } 70 71 // Summary output 72 if !quiet { 73 println!("Selected options: [command line option] : <current value>"); 74 println!(" Ceed specification [-c] : {}", ceed_spec); 75 println!(" Mesh dimension [-d] : {}", dim); 76 println!(" Mesh degree [-m] : {}", mesh_degree); 77 println!(" Solution degree [-p] : {}", solution_degree); 78 println!(" Num. 1D quadr. pts [-q] : {}", num_qpts); 79 println!(" Approx. # unknowns [-s] : {}", problem_size); 80 println!( 81 " QFunction source [-g] : {}", 82 if gallery { "gallery" } else { "user closure" } 83 ); 84 } 85 86 // Initalize ceed context 87 let ceed = Ceed::init(&ceed_spec); 88 89 // Mesh and solution bases 90 let basis_mesh = 91 ceed.basis_tensor_H1_Lagrange(dim, ncomp_x, mesh_degree + 1, num_qpts, QuadMode::Gauss)?; 92 let basis_solution = 93 ceed.basis_tensor_H1_Lagrange(dim, 1, solution_degree + 1, num_qpts, QuadMode::Gauss)?; 94 95 // Determine mesh size from approximate problem size 96 let num_xyz = mesh::cartesian_mesh_size(dim, solution_degree, problem_size); 97 if !quiet { 98 print!("\nMesh size : nx = {}", num_xyz[0]); 99 if dim > 1 { 100 print!(", ny = {}", num_xyz[1]); 101 } 102 if dim > 2 { 103 print!(", nz = {}", num_xyz[2]); 104 } 105 print!("\n"); 106 } 107 108 // Build ElemRestriction objects describing the mesh and solution discrete 109 // representations 110 let (restr_mesh, _) = 111 mesh::build_cartesian_restriction(&ceed, dim, num_xyz, mesh_degree, ncomp_x, num_qpts)?; 112 let (restr_solution, restr_qdata) = 113 mesh::build_cartesian_restriction(&ceed, dim, num_xyz, solution_degree, 1, num_qpts)?; 114 let mesh_size = restr_mesh.lvector_size(); 115 let solution_size = restr_solution.lvector_size(); 116 if !quiet { 117 println!("Number of mesh nodes : {}", mesh_size / dim); 118 println!("Number of solution nodes : {}", solution_size); 119 } 120 121 // Create a Vector with the mesh coordinates 122 let mut mesh_coords = mesh::cartesian_mesh_coords(&ceed, dim, num_xyz, mesh_degree, mesh_size)?; 123 124 // Apply a transformation to the mesh coordinates 125 let exact_volume = transform::transform_mesh_coordinates(dim, mesh_size, &mut mesh_coords)?; 126 127 // QFunction that builds the quadrature data for the mass operator 128 // -- QFunction from user closure 129 let build_mass = move |[jacobian, weights, ..]: QFunctionInputs, 130 [qdata, ..]: QFunctionOutputs| { 131 // Build quadrature data 132 match dim { 133 1 => qdata 134 .iter_mut() 135 .zip(jacobian.iter().zip(weights.iter())) 136 .for_each(|(qdata, (j, weight))| *qdata = j * weight), 137 2 => { 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] * jacobian[i + q * 3] 142 - jacobian[i + q * 1] * jacobian[i + q * 2]) 143 * weight 144 }, 145 ); 146 } 147 3 => { 148 let q = qdata.len(); 149 qdata.iter_mut().zip(weights.iter()).enumerate().for_each( 150 |(i, (qdata, weight))| { 151 *qdata = (jacobian[i + q * 0] 152 * (jacobian[i + q * 4] * jacobian[i + q * 8] 153 - jacobian[i + q * 5] * jacobian[i + q * 7]) 154 - jacobian[i + q * 1] 155 * (jacobian[i + q * 3] * jacobian[i + q * 8] 156 - jacobian[i + q * 5] * jacobian[i + q * 6]) 157 + jacobian[i + q * 2] 158 * (jacobian[i + q * 3] * jacobian[i + q * 7] 159 - jacobian[i + q * 4] * jacobian[i + q * 6])) 160 * weight 161 }, 162 ); 163 } 164 _ => unreachable!(), 165 }; 166 167 // Return clean error code 168 0 169 }; 170 let qf_build_closure = ceed 171 .q_function_interior(1, Box::new(build_mass))? 172 .input("dx", ncomp_x * dim, EvalMode::Grad)? 173 .input("weights", 1, EvalMode::Weight)? 174 .output("qdata", 1, EvalMode::None)?; 175 // -- QFunction from gallery 176 let qf_build_named = { 177 let name = format!("Mass{}DBuild", dim); 178 ceed.q_function_interior_by_name(&name)? 179 }; 180 // -- QFunction for use with Operator 181 let qf_build = if gallery { 182 QFunctionOpt::SomeQFunctionByName(&qf_build_named) 183 } else { 184 QFunctionOpt::SomeQFunction(&qf_build_closure) 185 }; 186 187 // Operator that build the quadrature data for the mass operator 188 let op_build = ceed 189 .operator(qf_build, QFunctionOpt::None, QFunctionOpt::None)? 190 .field("dx", &restr_mesh, &basis_mesh, VectorOpt::Active)? 191 .field( 192 "weights", 193 ElemRestrictionOpt::None, 194 &basis_mesh, 195 VectorOpt::None, 196 )? 197 .field( 198 "qdata", 199 &restr_qdata, 200 BasisOpt::Collocated, 201 VectorOpt::Active, 202 )?; 203 204 // Compute the quadrature data for the mass operator 205 let elem_qpts = num_qpts.pow(dim as u32); 206 let num_elem: usize = num_xyz.iter().take(dim).product(); 207 let mut qdata = ceed.vector(num_elem * elem_qpts)?; 208 op_build.apply(&mesh_coords, &mut qdata)?; 209 210 // QFunction that applies the mass operator 211 // -- QFunction from user closure 212 let apply_mass = |[u, qdata, ..]: QFunctionInputs, [v, ..]: QFunctionOutputs| { 213 // Apply mass operator 214 v.iter_mut() 215 .zip(u.iter().zip(qdata.iter())) 216 .for_each(|(v, (u, w))| *v = u * w); 217 // Return clean error code 218 0 219 }; 220 let qf_mass_closure = ceed 221 .q_function_interior(1, Box::new(apply_mass))? 222 .input("u", 1, EvalMode::Interp)? 223 .input("qdata", 1, EvalMode::None)? 224 .output("v", 1, EvalMode::Interp)?; 225 // -- QFunction from gallery 226 let qf_mass_named = ceed.q_function_interior_by_name("MassApply")?; 227 // -- QFunction for use with Operator 228 let qf_mass = if gallery { 229 QFunctionOpt::SomeQFunctionByName(&qf_mass_named) 230 } else { 231 QFunctionOpt::SomeQFunction(&qf_mass_closure) 232 }; 233 234 // Mass Operator 235 let op_mass = ceed 236 .operator(qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 237 .field("u", &restr_solution, &basis_solution, VectorOpt::Active)? 238 .field("qdata", &restr_qdata, BasisOpt::Collocated, &qdata)? 239 .field("v", &restr_solution, &basis_solution, VectorOpt::Active)?; 240 241 // Solution vectors 242 let u = ceed.vector_from_slice(&vec![1.0; solution_size])?; 243 let mut v = ceed.vector(solution_size)?; 244 245 // Apply the mass operator 246 op_mass.apply(&u, &mut v)?; 247 248 // Compute the mesh volume 249 let volume: Scalar = v.view()?.iter().sum(); 250 251 // Output results 252 if !quiet { 253 println!("Exact mesh volume : {:.12}", exact_volume); 254 println!("Computed mesh volume : {:.12}", volume); 255 println!( 256 "Volume error : {:.12e}", 257 volume - exact_volume 258 ); 259 } 260 let tolerance = match dim { 261 1 => 100.0 * libceed::EPSILON, 262 _ => 1E-5, 263 }; 264 let error = (volume - exact_volume).abs(); 265 if error > tolerance { 266 println!("Volume error too large: {:.12e}", error); 267 return Err(libceed::Error { 268 message: format!( 269 "Volume error too large - expected: {:.12e}, actual: {:.12e}", 270 tolerance, error 271 ), 272 }); 273 } 274 Ok(()) 275 } 276 277 // ---------------------------------------------------------------------------- 278 // Tests 279 // ---------------------------------------------------------------------------- 280 #[cfg(test)] 281 mod tests { 282 use super::*; 283 284 #[test] 285 fn example_1_1d() { 286 let options = opt::Opt { 287 ceed_spec: "/cpu/self/ref/serial".to_string(), 288 dim: 1, 289 mesh_degree: 4, 290 solution_degree: 4, 291 num_qpts: 6, 292 problem_size_requested: -1, 293 test: true, 294 quiet: true, 295 gallery: false, 296 }; 297 assert!(example_1(options).is_ok()); 298 } 299 300 #[test] 301 fn example_1_2d() { 302 let options = opt::Opt { 303 ceed_spec: "/cpu/self/ref/serial".to_string(), 304 dim: 2, 305 mesh_degree: 4, 306 solution_degree: 4, 307 num_qpts: 6, 308 problem_size_requested: -1, 309 test: true, 310 quiet: true, 311 gallery: false, 312 }; 313 assert!(example_1(options).is_ok()); 314 } 315 316 #[test] 317 fn example_1_3d() { 318 let options = opt::Opt { 319 ceed_spec: "/cpu/self/ref/serial".to_string(), 320 dim: 3, 321 mesh_degree: 4, 322 solution_degree: 4, 323 num_qpts: 6, 324 problem_size_requested: -1, 325 test: true, 326 quiet: false, 327 gallery: false, 328 }; 329 assert!(example_1(options).is_ok()); 330 } 331 332 #[test] 333 fn example_1_1d_gallery() { 334 let options = opt::Opt { 335 ceed_spec: "/cpu/self/ref/serial".to_string(), 336 dim: 1, 337 mesh_degree: 4, 338 solution_degree: 4, 339 num_qpts: 6, 340 problem_size_requested: -1, 341 test: true, 342 quiet: true, 343 gallery: true, 344 }; 345 assert!(example_1(options).is_ok()); 346 } 347 348 #[test] 349 fn example_1_2d_gallery() { 350 let options = opt::Opt { 351 ceed_spec: "/cpu/self/ref/serial".to_string(), 352 dim: 2, 353 mesh_degree: 4, 354 solution_degree: 4, 355 num_qpts: 6, 356 problem_size_requested: -1, 357 test: true, 358 quiet: true, 359 gallery: true, 360 }; 361 assert!(example_1(options).is_ok()); 362 } 363 364 #[test] 365 fn example_1_3d_gallery() { 366 let options = opt::Opt { 367 ceed_spec: "/cpu/self/ref/serial".to_string(), 368 dim: 3, 369 mesh_degree: 4, 370 solution_degree: 4, 371 num_qpts: 6, 372 problem_size_requested: -1, 373 test: true, 374 quiet: true, 375 gallery: true, 376 }; 377 assert!(example_1(options).is_ok()); 378 } 379 } 380 381 // ---------------------------------------------------------------------------- 382