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