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 .check()?; 204 205 // Compute the quadrature data for the mass operator 206 let elem_qpts = num_qpts.pow(dim as u32); 207 let num_elem: usize = num_xyz.iter().take(dim).product(); 208 let mut qdata = ceed.vector(num_elem * elem_qpts)?; 209 op_build.apply(&mesh_coords, &mut qdata)?; 210 211 // QFunction that applies the mass operator 212 // -- QFunction from user closure 213 let apply_mass = |[u, qdata, ..]: QFunctionInputs, [v, ..]: QFunctionOutputs| { 214 // Apply mass operator 215 v.iter_mut() 216 .zip(u.iter().zip(qdata.iter())) 217 .for_each(|(v, (u, w))| *v = u * w); 218 // Return clean error code 219 0 220 }; 221 let qf_mass_closure = ceed 222 .q_function_interior(1, Box::new(apply_mass))? 223 .input("u", 1, EvalMode::Interp)? 224 .input("qdata", 1, EvalMode::None)? 225 .output("v", 1, EvalMode::Interp)?; 226 // -- QFunction from gallery 227 let qf_mass_named = ceed.q_function_interior_by_name("MassApply")?; 228 // -- QFunction for use with Operator 229 let qf_mass = if gallery { 230 QFunctionOpt::SomeQFunctionByName(&qf_mass_named) 231 } else { 232 QFunctionOpt::SomeQFunction(&qf_mass_closure) 233 }; 234 235 // Mass Operator 236 let op_mass = ceed 237 .operator(qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 238 .field("u", &restr_solution, &basis_solution, VectorOpt::Active)? 239 .field("qdata", &restr_qdata, BasisOpt::Collocated, &qdata)? 240 .field("v", &restr_solution, &basis_solution, VectorOpt::Active)? 241 .check()?; 242 243 // Solution vectors 244 let u = ceed.vector_from_slice(&vec![1.0; solution_size])?; 245 let mut v = ceed.vector(solution_size)?; 246 247 // Apply the mass operator 248 op_mass.apply(&u, &mut v)?; 249 250 // Compute the mesh volume 251 let volume: Scalar = v.view()?.iter().sum(); 252 253 // Output results 254 if !quiet { 255 println!("Exact mesh volume : {:.12}", exact_volume); 256 println!("Computed mesh volume : {:.12}", volume); 257 println!( 258 "Volume error : {:.12e}", 259 volume - exact_volume 260 ); 261 } 262 let tolerance = match dim { 263 1 => 100.0 * libceed::EPSILON, 264 _ => 1E-5, 265 }; 266 let error = (volume - exact_volume).abs(); 267 if error > tolerance { 268 println!("Volume error too large: {:.12e}", error); 269 return Err(libceed::Error { 270 message: format!( 271 "Volume error too large - expected: {:.12e}, actual: {:.12e}", 272 tolerance, error 273 ), 274 }); 275 } 276 Ok(()) 277 } 278 279 // ---------------------------------------------------------------------------- 280 // Tests 281 // ---------------------------------------------------------------------------- 282 #[cfg(test)] 283 mod tests { 284 use super::*; 285 286 #[test] 287 fn example_1_1d() { 288 let options = opt::Opt { 289 ceed_spec: "/cpu/self/ref/serial".to_string(), 290 dim: 1, 291 mesh_degree: 4, 292 solution_degree: 4, 293 num_qpts: 6, 294 problem_size_requested: -1, 295 test: true, 296 quiet: true, 297 gallery: false, 298 }; 299 assert!(example_1(options).is_ok()); 300 } 301 302 #[test] 303 fn example_1_2d() { 304 let options = opt::Opt { 305 ceed_spec: "/cpu/self/ref/serial".to_string(), 306 dim: 2, 307 mesh_degree: 4, 308 solution_degree: 4, 309 num_qpts: 6, 310 problem_size_requested: -1, 311 test: true, 312 quiet: true, 313 gallery: false, 314 }; 315 assert!(example_1(options).is_ok()); 316 } 317 318 #[test] 319 fn example_1_3d() { 320 let options = opt::Opt { 321 ceed_spec: "/cpu/self/ref/serial".to_string(), 322 dim: 3, 323 mesh_degree: 4, 324 solution_degree: 4, 325 num_qpts: 6, 326 problem_size_requested: -1, 327 test: true, 328 quiet: false, 329 gallery: false, 330 }; 331 assert!(example_1(options).is_ok()); 332 } 333 334 #[test] 335 fn example_1_1d_gallery() { 336 let options = opt::Opt { 337 ceed_spec: "/cpu/self/ref/serial".to_string(), 338 dim: 1, 339 mesh_degree: 4, 340 solution_degree: 4, 341 num_qpts: 6, 342 problem_size_requested: -1, 343 test: true, 344 quiet: true, 345 gallery: true, 346 }; 347 assert!(example_1(options).is_ok()); 348 } 349 350 #[test] 351 fn example_1_2d_gallery() { 352 let options = opt::Opt { 353 ceed_spec: "/cpu/self/ref/serial".to_string(), 354 dim: 2, 355 mesh_degree: 4, 356 solution_degree: 4, 357 num_qpts: 6, 358 problem_size_requested: -1, 359 test: true, 360 quiet: true, 361 gallery: true, 362 }; 363 assert!(example_1(options).is_ok()); 364 } 365 366 #[test] 367 fn example_1_3d_gallery() { 368 let options = opt::Opt { 369 ceed_spec: "/cpu/self/ref/serial".to_string(), 370 dim: 3, 371 mesh_degree: 4, 372 solution_degree: 4, 373 num_qpts: 6, 374 problem_size_requested: -1, 375 test: true, 376 quiet: true, 377 gallery: true, 378 }; 379 assert!(example_1(options).is_ok()); 380 } 381 } 382 383 // ---------------------------------------------------------------------------- 384