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