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