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 = ceed 91 .basis_tensor_H1_Lagrange(dim, ncomp_x, mesh_degree + 1, num_qpts, QuadMode::Gauss) 92 .unwrap(); 93 let basis_solution = ceed 94 .basis_tensor_H1_Lagrange(dim, 1, solution_degree + 1, num_qpts, QuadMode::Gauss) 95 .unwrap(); 96 97 // Determine mesh size from approximate problem size 98 let num_xyz = mesh::cartesian_mesh_size(dim, solution_degree, problem_size); 99 if !quiet { 100 print!("\nMesh size : nx = {}", num_xyz[0]); 101 if dim > 1 { 102 print!(", ny = {}", num_xyz[1]); 103 } 104 if dim > 2 { 105 print!(", nz = {}", num_xyz[2]); 106 } 107 print!("\n"); 108 } 109 110 // Build ElemRestriction objects describing the mesh and solution discrete 111 // representations 112 let (restr_mesh, _) = 113 mesh::build_cartesian_restriction(&ceed, dim, num_xyz, mesh_degree, ncomp_x, num_qpts)?; 114 let (restr_solution, restr_qdata) = 115 mesh::build_cartesian_restriction(&ceed, dim, num_xyz, solution_degree, 1, num_qpts)?; 116 let mesh_size = restr_mesh.lvector_size(); 117 let solution_size = restr_solution.lvector_size(); 118 if !quiet { 119 println!("Number of mesh nodes : {}", mesh_size / dim); 120 println!("Number of solution nodes : {}", solution_size); 121 } 122 123 // Create a Vector with the mesh coordinates 124 let mut mesh_coords = mesh::cartesian_mesh_coords(&ceed, dim, num_xyz, mesh_degree, mesh_size)?; 125 126 // Apply a transformation to the mesh coordinates 127 let exact_volume = transform::transform_mesh_coordinates(dim, mesh_size, &mut mesh_coords)?; 128 129 // QFunction that builds the quadrature data for the mass operator 130 // -- QFunction from user closure 131 let build_mass = move |[jacobian, weights, ..]: QFunctionInputs, 132 [qdata, ..]: QFunctionOutputs| { 133 // Build quadrature data 134 match dim { 135 1 => qdata 136 .iter_mut() 137 .zip(jacobian.iter().zip(weights.iter())) 138 .for_each(|(qdata, (j, weight))| *qdata = j * weight), 139 2 => { 140 let q = qdata.len(); 141 qdata.iter_mut().zip(weights.iter()).enumerate().for_each( 142 |(i, (qdata, weight))| { 143 *qdata = (jacobian[i + q * 0] * jacobian[i + q * 3] 144 - jacobian[i + q * 1] * jacobian[i + q * 2]) 145 * weight 146 }, 147 ); 148 } 149 3 => { 150 let q = qdata.len(); 151 qdata.iter_mut().zip(weights.iter()).enumerate().for_each( 152 |(i, (qdata, weight))| { 153 *qdata = (jacobian[i + q * 0] 154 * (jacobian[i + q * 4] * jacobian[i + q * 8] 155 - jacobian[i + q * 5] * jacobian[i + q * 7]) 156 - jacobian[i + q * 1] 157 * (jacobian[i + q * 3] * jacobian[i + q * 8] 158 - jacobian[i + q * 5] * jacobian[i + q * 6]) 159 + jacobian[i + q * 2] 160 * (jacobian[i + q * 3] * jacobian[i + q * 7] 161 - jacobian[i + q * 4] * jacobian[i + q * 6])) 162 * weight 163 }, 164 ); 165 } 166 _ => unreachable!(), 167 }; 168 169 // Return clean error code 170 0 171 }; 172 let qf_build_closure = ceed 173 .q_function_interior(1, Box::new(build_mass))? 174 .input("dx", ncomp_x * dim, EvalMode::Grad)? 175 .input("weights", 1, EvalMode::Weight)? 176 .output("qdata", 1, EvalMode::None)?; 177 // -- QFunction from gallery 178 let qf_build_named = { 179 let name = format!("Mass{}DBuild", dim); 180 ceed.q_function_interior_by_name(&name)? 181 }; 182 // -- QFunction for use with Operator 183 let qf_build = if gallery { 184 QFunctionOpt::SomeQFunctionByName(&qf_build_named) 185 } else { 186 QFunctionOpt::SomeQFunction(&qf_build_closure) 187 }; 188 189 // Operator that build the quadrature data for the mass operator 190 let op_build = ceed 191 .operator(qf_build, QFunctionOpt::None, QFunctionOpt::None)? 192 .field("dx", &restr_mesh, &basis_mesh, VectorOpt::Active)? 193 .field( 194 "weights", 195 ElemRestrictionOpt::None, 196 &basis_mesh, 197 VectorOpt::None, 198 )? 199 .field( 200 "qdata", 201 &restr_qdata, 202 BasisOpt::Collocated, 203 VectorOpt::Active, 204 )?; 205 206 // Compute the quadrature data for the mass operator 207 let elem_qpts = num_qpts.pow(dim as u32); 208 let num_elem: usize = num_xyz.iter().take(dim).product(); 209 let mut qdata = ceed.vector(num_elem * elem_qpts)?; 210 op_build.apply(&mesh_coords, &mut qdata)?; 211 212 // QFunction that applies the mass operator 213 // -- QFunction from user closure 214 let apply_mass = |[u, qdata, ..]: QFunctionInputs, [v, ..]: QFunctionOutputs| { 215 // Apply mass operator 216 v.iter_mut() 217 .zip(u.iter().zip(qdata.iter())) 218 .for_each(|(v, (u, w))| *v = u * w); 219 // Return clean error code 220 0 221 }; 222 let qf_mass_closure = ceed 223 .q_function_interior(1, Box::new(apply_mass))? 224 .input("u", 1, EvalMode::Interp)? 225 .input("qdata", 1, EvalMode::None)? 226 .output("v", 1, EvalMode::Interp)?; 227 // -- QFunction from gallery 228 let qf_mass_named = ceed.q_function_interior_by_name("MassApply")?; 229 // -- QFunction for use with Operator 230 let qf_mass = if gallery { 231 QFunctionOpt::SomeQFunctionByName(&qf_mass_named) 232 } else { 233 QFunctionOpt::SomeQFunction(&qf_mass_closure) 234 }; 235 236 // Mass Operator 237 let op_mass = ceed 238 .operator(qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 239 .field("u", &restr_solution, &basis_solution, VectorOpt::Active)? 240 .field("qdata", &restr_qdata, BasisOpt::Collocated, &qdata)? 241 .field("v", &restr_solution, &basis_solution, VectorOpt::Active)?; 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