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 2 18 // 19 // This example illustrates a simple usage of libCEED to compute the surface 20 // area of a 3D body using matrix-free application of a diffusion operator. 21 // Arbitrary mesh and solution degrees in 1D, 2D and 3D are supported from the 22 // same code. 23 // 24 // The example has no dependencies, and is designed to be self-contained. For 25 // additional examples that use external discretization libraries (MFEM, PETSc, 26 // etc.) see the subdirectories in libceed/examples. 27 // 28 // All libCEED objects use a Ceed device object constructed based on a command 29 // line argument (-ceed). 30 31 use libceed::{prelude::*, Ceed}; 32 use mesh; 33 use structopt::StructOpt; 34 35 mod opt; 36 mod transform; 37 38 // ---------------------------------------------------------------------------- 39 // Example 2 40 // ---------------------------------------------------------------------------- 41 #[cfg(not(tarpaulin_include))] 42 fn main() -> Result<(), libceed::CeedError> { 43 let options = opt::Opt::from_args(); 44 example_2(options) 45 } 46 47 fn example_2(options: opt::Opt) -> Result<(), libceed::CeedError> { 48 // Process command line arguments 49 let opt::Opt { 50 ceed_spec, 51 dim, 52 mesh_degree, 53 solution_degree, 54 num_qpts, 55 problem_size_requested, 56 test, 57 quiet, 58 gallery, 59 } = options; 60 assert!(dim >= 1 && dim <= 3); 61 assert!(mesh_degree >= 1); 62 assert!(solution_degree >= 1); 63 assert!(num_qpts >= 1); 64 let ncomp_x = dim; 65 let problem_size: i64; 66 if problem_size_requested < 0 { 67 problem_size = if test { 68 16 * 16 * (dim * dim) as i64 69 } else { 70 256 * 1024 71 }; 72 } else { 73 problem_size = problem_size_requested; 74 } 75 76 // Summary output 77 if !quiet { 78 println!("Selected options: [command line option] : <current value>"); 79 println!(" Ceed specification [-c] : {}", ceed_spec); 80 println!(" Mesh dimension [-d] : {}", dim); 81 println!(" Mesh degree [-m] : {}", mesh_degree); 82 println!(" Solution degree [-p] : {}", solution_degree); 83 println!(" Num. 1D quadr. pts [-q] : {}", num_qpts); 84 println!(" Approx. # unknowns [-s] : {}", problem_size); 85 println!( 86 " QFunction source [-g] : {}", 87 if gallery { "gallery" } else { "user closure" } 88 ); 89 } 90 91 // Initalize ceed context 92 let ceed = Ceed::init(&ceed_spec); 93 94 // Mesh and solution bases 95 let basis_mesh = ceed 96 .basis_tensor_H1_Lagrange(dim, ncomp_x, mesh_degree + 1, num_qpts, QuadMode::Gauss) 97 .unwrap(); 98 let basis_solution = ceed 99 .basis_tensor_H1_Lagrange(dim, 1, solution_degree + 1, num_qpts, QuadMode::Gauss) 100 .unwrap(); 101 102 // Determine mesh size from approximate problem size 103 let num_xyz = mesh::cartesian_mesh_size(dim, solution_degree, problem_size); 104 if !quiet { 105 print!("\nMesh size : nx = {}", num_xyz[0]); 106 if dim > 1 { 107 print!(", ny = {}", num_xyz[1]); 108 } 109 if dim > 2 { 110 print!(", nz = {}", num_xyz[2]); 111 } 112 print!("\n"); 113 } 114 115 // Build ElemRestriction objects describing the mesh and solution discrete 116 // representations 117 let (restr_mesh, _) = 118 mesh::build_cartesian_restriction(&ceed, dim, num_xyz, mesh_degree, ncomp_x, num_qpts)?; 119 let (_, restr_qdata) = mesh::build_cartesian_restriction( 120 &ceed, 121 dim, 122 num_xyz, 123 solution_degree, 124 dim * (dim + 1) / 2, 125 num_qpts, 126 )?; 127 128 let (restr_solution, _) = 129 mesh::build_cartesian_restriction(&ceed, dim, num_xyz, solution_degree, 1, num_qpts)?; 130 let mesh_size = restr_mesh.lvector_size(); 131 let solution_size = restr_solution.lvector_size(); 132 if !quiet { 133 println!("Number of mesh nodes : {}", mesh_size / dim); 134 println!("Number of solution nodes : {}", solution_size); 135 } 136 137 // Create a Vector with the mesh coordinates 138 let mut mesh_coords = mesh::cartesian_mesh_coords(&ceed, dim, num_xyz, mesh_degree, mesh_size)?; 139 140 // Apply a transformation to the mesh coordinates 141 let exact_area = transform::transform_mesh_coordinates(dim, &mut mesh_coords); 142 143 // QFunction that builds the quadrature data for the diff operator 144 // -- QFunction from user closure 145 let build_diff = move |[jacobian, weights, ..]: QFunctionInputs, 146 [qdata, ..]: QFunctionOutputs| { 147 // Build quadrature data 148 match dim { 149 1 => qdata 150 .iter_mut() 151 .zip(jacobian.iter().zip(weights.iter())) 152 .for_each(|(qdata, (j, weight))| *qdata = weight / j), 153 2 => { 154 let q = qdata.len() / 3; 155 for i in 0..q { 156 let j11 = jacobian[i + q * 0]; 157 let j21 = jacobian[i + q * 1]; 158 let j12 = jacobian[i + q * 2]; 159 let j22 = jacobian[i + q * 3]; 160 let qw = weights[i] / (j11 * j22 - j21 * j12); 161 qdata[i + q * 0] = qw * (j12 * j12 + j22 * j22); 162 qdata[i + q * 1] = qw * (j11 * j11 + j21 * j21); 163 qdata[i + q * 2] = -qw * (j11 * j12 + j21 * j22); 164 } 165 } 166 3 => { 167 let q = qdata.len() / 6; 168 for i in 0..q { 169 let mut a = [0.0; 9]; 170 for j in 0..3 { 171 for k in 0..3 { 172 a[k * 3 + j] = jacobian[i + q * ((j + 1) % 3 + 3 * ((k + 1) % 3))] 173 * jacobian[i + q * ((j + 2) % 3 + 3 * ((k + 2) % 3))] 174 - jacobian[i + q * ((j + 1) % 3 + 3 * ((k + 2) % 3))] 175 * jacobian[i + q * ((j + 2) % 3 + 3 * ((k + 1) % 3))]; 176 } 177 } 178 let qw = weights[i] 179 / (jacobian[i + q * 0] * a[0 * 3 + 0] 180 + jacobian[i + q * 1] * a[1 * 3 + 1] 181 + jacobian[i + q * 2] * a[2 * 3 + 2]); 182 qdata[i + q * 0] = qw 183 * (a[0 * 3 + 0] * a[0 * 3 + 0] 184 + a[0 * 3 + 1] * a[0 * 3 + 1] 185 + a[0 * 3 + 2] * a[0 * 3 + 2]); 186 qdata[i + q * 1] = qw 187 * (a[1 * 3 + 0] * a[1 * 3 + 0] 188 + a[1 * 3 + 1] * a[1 * 3 + 1] 189 + a[1 * 3 + 2] * a[1 * 3 + 2]); 190 qdata[i + q * 2] = qw 191 * (a[2 * 3 + 0] * a[2 * 3 + 0] 192 + a[2 * 3 + 1] * a[2 * 3 + 1] 193 + a[2 * 3 + 2] * a[2 * 3 + 2]); 194 qdata[i + q * 3] = qw 195 * (a[1 * 3 + 0] * a[2 * 3 + 0] 196 + a[1 * 3 + 1] * a[2 * 3 + 1] 197 + a[1 * 3 + 2] * a[2 * 3 + 2]); 198 qdata[i + q * 4] = qw 199 * (a[0 * 3 + 0] * a[2 * 3 + 0] 200 + a[0 * 3 + 1] * a[2 * 3 + 1] 201 + a[0 * 3 + 2] * a[2 * 3 + 2]); 202 qdata[i + q * 5] = qw 203 * (a[0 * 3 + 0] * a[1 * 3 + 0] 204 + a[0 * 3 + 1] * a[1 * 3 + 1] 205 + a[0 * 3 + 2] * a[1 * 3 + 2]); 206 } 207 } 208 _ => unreachable!(), 209 }; 210 211 // Return clean error code 212 0 213 }; 214 let qf_build_closure = ceed 215 .q_function_interior(1, Box::new(build_diff))? 216 .input("dx", ncomp_x * dim, EvalMode::Grad)? 217 .input("weights", 1, EvalMode::Weight)? 218 .output("qdata", dim * (dim + 1) / 2, EvalMode::None)?; 219 // -- QFunction from gallery 220 let qf_build_named = { 221 let name = format!("Poisson{}DBuild", dim); 222 ceed.q_function_interior_by_name(&name)? 223 }; 224 // -- QFunction for use with Operator 225 let qf_build = if gallery { 226 QFunctionOpt::SomeQFunctionByName(&qf_build_named) 227 } else { 228 QFunctionOpt::SomeQFunction(&qf_build_closure) 229 }; 230 231 // Operator that build the quadrature data for the diff operator 232 let op_build = ceed 233 .operator(qf_build, QFunctionOpt::None, QFunctionOpt::None)? 234 .field("dx", &restr_mesh, &basis_mesh, VectorOpt::Active)? 235 .field( 236 "weights", 237 ElemRestrictionOpt::None, 238 &basis_mesh, 239 VectorOpt::None, 240 )? 241 .field( 242 "qdata", 243 &restr_qdata, 244 BasisOpt::Collocated, 245 VectorOpt::Active, 246 )?; 247 248 // Compute the quadrature data for the diff operator 249 let elem_qpts = num_qpts.pow(dim as u32); 250 let num_elem: usize = num_xyz.iter().take(dim).product(); 251 let mut qdata = ceed.vector(num_elem * elem_qpts * dim * (dim + 1) / 2)?; 252 op_build.apply(&mesh_coords, &mut qdata)?; 253 254 // QFunction that applies the diff operator 255 // -- QFunction from user closure 256 let apply_diff = move |[ug, qdata, ..]: QFunctionInputs, [vg, ..]: QFunctionOutputs| { 257 // Apply diffusion operator 258 match dim { 259 1 => vg 260 .iter_mut() 261 .zip(ug.iter().zip(qdata.iter())) 262 .for_each(|(vg, (ug, w))| *vg = ug * w), 263 2 => { 264 let q = qdata.len() / 3; 265 for i in 0..q { 266 let du = [ug[i + q * 0], ug[i + q * 1]]; 267 let dxdxdxdx_t = [ 268 [qdata[i + 0 * q], qdata[i + 2 * q]], 269 [qdata[i + 2 * q], qdata[i + 1 * q]], 270 ]; 271 for j in 0..2 { 272 vg[i + j * q] = du[0] * dxdxdxdx_t[0][j] + du[1] * dxdxdxdx_t[1][j]; 273 } 274 } 275 } 276 3 => { 277 let q = qdata.len() / 6; 278 for i in 0..q { 279 let du = [ug[i + q * 0], ug[i + q * 1], ug[i + q * 2]]; 280 let dxdxdxdx_t = [ 281 [qdata[i + 0 * q], qdata[i + 5 * q], qdata[i + 4 * q]], 282 [qdata[i + 5 * q], qdata[i + 1 * q], qdata[i + 3 * q]], 283 [qdata[i + 4 * q], qdata[i + 3 * q], qdata[i + 2 * q]], 284 ]; 285 for j in 0..3 { 286 vg[i + j * q] = du[0] * dxdxdxdx_t[0][j] 287 + du[1] * dxdxdxdx_t[1][j] 288 + du[2] * dxdxdxdx_t[2][j]; 289 } 290 } 291 } 292 _ => unreachable!(), 293 }; 294 295 // Return clean error code 296 0 297 }; 298 let qf_diff_closure = ceed 299 .q_function_interior(1, Box::new(apply_diff))? 300 .input("du", dim, EvalMode::Grad)? 301 .input("qdata", dim * (dim + 1) / 2, EvalMode::None)? 302 .output("dv", dim, EvalMode::Grad)?; 303 // -- QFunction from gallery 304 let qf_diff_named = { 305 let name = format!("Poisson{}DApply", dim); 306 ceed.q_function_interior_by_name(&name)? 307 }; 308 // -- QFunction for use with Operator 309 let qf_diff = if gallery { 310 QFunctionOpt::SomeQFunctionByName(&qf_diff_named) 311 } else { 312 QFunctionOpt::SomeQFunction(&qf_diff_closure) 313 }; 314 315 // Diff Operator 316 let op_diff = ceed 317 .operator(qf_diff, QFunctionOpt::None, QFunctionOpt::None)? 318 .field("du", &restr_solution, &basis_solution, VectorOpt::Active)? 319 .field("qdata", &restr_qdata, BasisOpt::Collocated, &qdata)? 320 .field("dv", &restr_solution, &basis_solution, VectorOpt::Active)?; 321 322 // Solution vectors 323 let mut u = ceed.vector(solution_size)?; 324 let mut v = ceed.vector(solution_size)?; 325 326 // Initialize u with sum of node coordinates 327 let coords = mesh_coords.view(); 328 u.view_mut().iter_mut().enumerate().for_each(|(i, u)| { 329 *u = (0..dim).map(|d| coords[i + d * solution_size]).sum(); 330 }); 331 332 // Apply the diff operator 333 op_diff.apply(&u, &mut v)?; 334 335 // Compute the mesh surface area 336 let area: f64 = v.view().iter().map(|v| (*v).abs()).sum(); 337 338 // Output results 339 if !quiet { 340 println!("Exact mesh surface area : {:.12}", exact_area); 341 println!("Computed mesh surface_area : {:.12}", area); 342 println!("Surface area error : {:.12e}", area - exact_area); 343 } 344 let tolerance = match dim { 345 1 => 1E-12, 346 _ => 1E-1, 347 }; 348 let error = (area - exact_area).abs(); 349 if error > tolerance { 350 println!("Volume error too large: {:.12e}", error); 351 return Err(libceed::CeedError { 352 message: format!( 353 "Volume error too large - expected: {:.12e}, actual: {:.12e}", 354 tolerance, error 355 ), 356 }); 357 } 358 Ok(()) 359 } 360 361 // ---------------------------------------------------------------------------- 362 // Tests 363 // ---------------------------------------------------------------------------- 364 #[cfg(test)] 365 mod tests { 366 use super::*; 367 368 #[test] 369 fn example_2_1d() { 370 let options = opt::Opt { 371 ceed_spec: "/cpu/self/ref/serial".to_string(), 372 dim: 1, 373 mesh_degree: 4, 374 solution_degree: 4, 375 num_qpts: 6, 376 problem_size_requested: -1, 377 test: true, 378 quiet: true, 379 gallery: false, 380 }; 381 assert!(example_2(options).is_ok()); 382 } 383 384 #[test] 385 fn example_2_2d() { 386 let options = opt::Opt { 387 ceed_spec: "/cpu/self/ref/serial".to_string(), 388 dim: 2, 389 mesh_degree: 4, 390 solution_degree: 4, 391 num_qpts: 6, 392 problem_size_requested: -1, 393 test: true, 394 quiet: true, 395 gallery: false, 396 }; 397 assert!(example_2(options).is_ok()); 398 } 399 400 #[test] 401 fn example_2_3d() { 402 let options = opt::Opt { 403 ceed_spec: "/cpu/self/ref/serial".to_string(), 404 dim: 3, 405 mesh_degree: 4, 406 solution_degree: 4, 407 num_qpts: 6, 408 problem_size_requested: -1, 409 test: true, 410 quiet: false, 411 gallery: false, 412 }; 413 assert!(example_2(options).is_ok()); 414 } 415 416 #[test] 417 fn example_2_1d_gallery() { 418 let options = opt::Opt { 419 ceed_spec: "/cpu/self/ref/serial".to_string(), 420 dim: 1, 421 mesh_degree: 4, 422 solution_degree: 4, 423 num_qpts: 6, 424 problem_size_requested: -1, 425 test: true, 426 quiet: true, 427 gallery: true, 428 }; 429 assert!(example_2(options).is_ok()); 430 } 431 432 #[test] 433 fn example_2_2d_gallery() { 434 let options = opt::Opt { 435 ceed_spec: "/cpu/self/ref/serial".to_string(), 436 dim: 2, 437 mesh_degree: 4, 438 solution_degree: 4, 439 num_qpts: 6, 440 problem_size_requested: -1, 441 test: true, 442 quiet: true, 443 gallery: true, 444 }; 445 assert!(example_2(options).is_ok()); 446 } 447 448 #[test] 449 fn example_2_3d_gallery() { 450 let options = opt::Opt { 451 ceed_spec: "/cpu/self/ref/serial".to_string(), 452 dim: 3, 453 mesh_degree: 4, 454 solution_degree: 4, 455 num_qpts: 6, 456 problem_size_requested: -1, 457 test: true, 458 quiet: true, 459 gallery: true, 460 }; 461 assert!(example_2(options).is_ok()); 462 } 463 } 464 465 // ---------------------------------------------------------------------------- 466