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<(), String> { 43 let options = opt::Opt::from_args(); 44 example_2(options) 45 } 46 47 fn example_2(options: opt::Opt) -> Result<(), String> { 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 = 96 ceed.basis_tensor_H1_Lagrange(dim, ncomp_x, mesh_degree + 1, num_qpts, QuadMode::Gauss); 97 let basis_solution = 98 ceed.basis_tensor_H1_Lagrange(dim, 1, solution_degree + 1, num_qpts, QuadMode::Gauss); 99 100 // Determine mesh size from approximate problem size 101 let num_xyz = mesh::cartesian_mesh_size(dim, solution_degree, problem_size); 102 if !quiet { 103 print!("\nMesh size : nx = {}", num_xyz[0]); 104 if dim > 1 { 105 print!(", ny = {}", num_xyz[1]); 106 } 107 if dim > 2 { 108 print!(", nz = {}", num_xyz[2]); 109 } 110 print!("\n"); 111 } 112 113 // Build ElemRestriction objects describing the mesh and solution discrete 114 // representations 115 let (restr_mesh, _) = 116 mesh::build_cartesian_restriction(&ceed, dim, num_xyz, mesh_degree, ncomp_x, num_qpts); 117 let (_, restr_qdata) = mesh::build_cartesian_restriction( 118 &ceed, 119 dim, 120 num_xyz, 121 solution_degree, 122 dim * (dim + 1) / 2, 123 num_qpts, 124 ); 125 126 let (restr_solution, _) = 127 mesh::build_cartesian_restriction(&ceed, dim, num_xyz, solution_degree, 1, num_qpts); 128 let mesh_size = restr_mesh.lvector_size(); 129 let solution_size = restr_solution.lvector_size(); 130 if !quiet { 131 println!("Number of mesh nodes : {}", mesh_size / dim); 132 println!("Number of solution nodes : {}", solution_size); 133 } 134 135 // Create a Vector with the mesh coordinates 136 let mut mesh_coords = mesh::cartesian_mesh_coords(&ceed, dim, num_xyz, mesh_degree, mesh_size); 137 138 // Apply a transformation to the mesh coordinates 139 let exact_area = transform::transform_mesh_coordinates(dim, &mut mesh_coords); 140 141 // QFunction that builds the quadrature data for the diff operator 142 // -- QFunction from user closure 143 let build_diff = move |[jacobian, weights, ..]: QFunctionInputs, 144 [qdata, ..]: QFunctionOutputs| { 145 // Build quadrature data 146 match dim { 147 1 => qdata 148 .iter_mut() 149 .zip(jacobian.iter().zip(weights.iter())) 150 .for_each(|(qdata, (j, weight))| *qdata = weight / j), 151 2 => { 152 let q = qdata.len() / 3; 153 for i in 0..q { 154 let j11 = jacobian[i + q * 0]; 155 let j21 = jacobian[i + q * 1]; 156 let j12 = jacobian[i + q * 2]; 157 let j22 = jacobian[i + q * 3]; 158 let qw = weights[i] / (j11 * j22 - j21 * j12); 159 qdata[i + q * 0] = qw * (j12 * j12 + j22 * j22); 160 qdata[i + q * 1] = qw * (j11 * j11 + j21 * j21); 161 qdata[i + q * 2] = -qw * (j11 * j12 + j21 * j22); 162 } 163 } 164 3 => { 165 let q = qdata.len() / 6; 166 for i in 0..q { 167 let mut a = [0.0; 9]; 168 for j in 0..3 { 169 for k in 0..3 { 170 a[k * 3 + j] = jacobian[i + q * ((j + 1) % 3 + 3 * ((k + 1) % 3))] 171 * jacobian[i + q * ((j + 2) % 3 + 3 * ((k + 2) % 3))] 172 - jacobian[i + q * ((j + 1) % 3 + 3 * ((k + 2) % 3))] 173 * jacobian[i + q * ((j + 2) % 3 + 3 * ((k + 1) % 3))]; 174 } 175 } 176 let qw = weights[i] 177 / (jacobian[i + q * 0] * a[0 * 3 + 0] 178 + jacobian[i + q * 1] * a[1 * 3 + 1] 179 + jacobian[i + q * 2] * a[2 * 3 + 2]); 180 qdata[i + q * 0] = qw 181 * (a[0 * 3 + 0] * a[0 * 3 + 0] 182 + a[0 * 3 + 1] * a[0 * 3 + 1] 183 + a[0 * 3 + 2] * a[0 * 3 + 2]); 184 qdata[i + q * 1] = qw 185 * (a[1 * 3 + 0] * a[1 * 3 + 0] 186 + a[1 * 3 + 1] * a[1 * 3 + 1] 187 + a[1 * 3 + 2] * a[1 * 3 + 2]); 188 qdata[i + q * 2] = qw 189 * (a[2 * 3 + 0] * a[2 * 3 + 0] 190 + a[2 * 3 + 1] * a[2 * 3 + 1] 191 + a[2 * 3 + 2] * a[2 * 3 + 2]); 192 qdata[i + q * 3] = qw 193 * (a[1 * 3 + 0] * a[2 * 3 + 0] 194 + a[1 * 3 + 1] * a[2 * 3 + 1] 195 + a[1 * 3 + 2] * a[2 * 3 + 2]); 196 qdata[i + q * 4] = qw 197 * (a[0 * 3 + 0] * a[2 * 3 + 0] 198 + a[0 * 3 + 1] * a[2 * 3 + 1] 199 + a[0 * 3 + 2] * a[2 * 3 + 2]); 200 qdata[i + q * 5] = qw 201 * (a[0 * 3 + 0] * a[1 * 3 + 0] 202 + a[0 * 3 + 1] * a[1 * 3 + 1] 203 + a[0 * 3 + 2] * a[1 * 3 + 2]); 204 } 205 } 206 _ => unreachable!(), 207 }; 208 209 // Return clean error code 210 0 211 }; 212 let qf_build_closure = ceed 213 .q_function_interior(1, Box::new(build_diff)) 214 .input("dx", ncomp_x * dim, EvalMode::Grad) 215 .input("weights", 1, EvalMode::Weight) 216 .output("qdata", dim * (dim + 1) / 2, EvalMode::None); 217 // -- QFunction from gallery 218 let qf_build_named = { 219 let name = format!("Poisson{}DBuild", dim); 220 ceed.q_function_interior_by_name(&name) 221 }; 222 // -- QFunction for use with Operator 223 let qf_build = if gallery { 224 QFunctionOpt::SomeQFunctionByName(&qf_build_named) 225 } else { 226 QFunctionOpt::SomeQFunction(&qf_build_closure) 227 }; 228 229 // Operator that build the quadrature data for the diff operator 230 let op_build = ceed 231 .operator(qf_build, QFunctionOpt::None, QFunctionOpt::None) 232 .field("dx", &restr_mesh, &basis_mesh, VectorOpt::Active) 233 .field( 234 "weights", 235 ElemRestrictionOpt::None, 236 &basis_mesh, 237 VectorOpt::None, 238 ) 239 .field( 240 "qdata", 241 &restr_qdata, 242 BasisOpt::Collocated, 243 VectorOpt::Active, 244 ); 245 246 // Compute the quadrature data for the diff operator 247 let elem_qpts = num_qpts.pow(dim as u32); 248 let num_elem: usize = num_xyz.iter().take(dim).product(); 249 let mut qdata = ceed.vector(num_elem * elem_qpts * dim * (dim + 1) / 2); 250 op_build.apply(&mesh_coords, &mut qdata); 251 252 // QFunction that applies the diff operator 253 // -- QFunction from user closure 254 let apply_diff = move |[ug, qdata, ..]: QFunctionInputs, [vg, ..]: QFunctionOutputs| { 255 // Apply diffusion operator 256 match dim { 257 1 => vg 258 .iter_mut() 259 .zip(ug.iter().zip(qdata.iter())) 260 .for_each(|(vg, (ug, w))| *vg = ug * w), 261 2 => { 262 let q = qdata.len() / 3; 263 for i in 0..q { 264 let du = [ug[i + q * 0], ug[i + q * 1]]; 265 let dxdxdxdx_t = [ 266 [qdata[i + 0 * q], qdata[i + 2 * q]], 267 [qdata[i + 2 * q], qdata[i + 1 * q]], 268 ]; 269 for j in 0..2 { 270 vg[i + j * q] = du[0] * dxdxdxdx_t[0][j] + du[1] * dxdxdxdx_t[1][j]; 271 } 272 } 273 } 274 3 => { 275 let q = qdata.len() / 6; 276 for i in 0..q { 277 let du = [ug[i + q * 0], ug[i + q * 1], ug[i + q * 2]]; 278 let dxdxdxdx_t = [ 279 [qdata[i + 0 * q], qdata[i + 5 * q], qdata[i + 4 * q]], 280 [qdata[i + 5 * q], qdata[i + 1 * q], qdata[i + 3 * q]], 281 [qdata[i + 4 * q], qdata[i + 3 * q], qdata[i + 2 * q]], 282 ]; 283 for j in 0..3 { 284 vg[i + j * q] = du[0] * dxdxdxdx_t[0][j] 285 + du[1] * dxdxdxdx_t[1][j] 286 + du[2] * dxdxdxdx_t[2][j]; 287 } 288 } 289 } 290 _ => unreachable!(), 291 }; 292 293 // Return clean error code 294 0 295 }; 296 let qf_diff_closure = ceed 297 .q_function_interior(1, Box::new(apply_diff)) 298 .input("du", dim, EvalMode::Grad) 299 .input("qdata", dim * (dim + 1) / 2, EvalMode::None) 300 .output("dv", dim, EvalMode::Grad); 301 // -- QFunction from gallery 302 let qf_diff_named = { 303 let name = format!("Poisson{}DApply", dim); 304 ceed.q_function_interior_by_name(&name) 305 }; 306 // -- QFunction for use with Operator 307 let qf_diff = if gallery { 308 QFunctionOpt::SomeQFunctionByName(&qf_diff_named) 309 } else { 310 QFunctionOpt::SomeQFunction(&qf_diff_closure) 311 }; 312 313 // Diff Operator 314 let op_diff = ceed 315 .operator(qf_diff, QFunctionOpt::None, QFunctionOpt::None) 316 .field("du", &restr_solution, &basis_solution, VectorOpt::Active) 317 .field("qdata", &restr_qdata, BasisOpt::Collocated, &qdata) 318 .field("dv", &restr_solution, &basis_solution, VectorOpt::Active); 319 320 // Solution vectors 321 let mut u = ceed.vector(solution_size); 322 let mut v = ceed.vector(solution_size); 323 324 // Initialize u with sum of node coordinates 325 let coords = mesh_coords.view(); 326 u.view_mut().iter_mut().enumerate().for_each(|(i, u)| { 327 *u = (0..dim).map(|d| coords[i + d * solution_size]).sum(); 328 }); 329 330 // Apply the diff operator 331 op_diff.apply(&u, &mut v); 332 333 // Compute the mesh surface area 334 let area: f64 = v.view().iter().map(|v| (*v).abs()).sum(); 335 336 // Output results 337 if !quiet { 338 println!("Exact mesh surface area : {:.12}", exact_area); 339 println!("Computed mesh surface_area : {:.12}", area); 340 println!("Surface area error : {:.12e}", area - exact_area); 341 } 342 let tolerance = match dim { 343 1 => 1E-12, 344 _ => 1E-1, 345 }; 346 let error = (area - exact_area).abs(); 347 if error > tolerance { 348 println!("Volume error too large: {:.12e}", error); 349 return Err(format!( 350 "Volume error too large - expected: {:.12e}, actual: {:.12e}", 351 tolerance, error 352 )); 353 } 354 Ok(()) 355 } 356 357 // ---------------------------------------------------------------------------- 358 // Tests 359 // ---------------------------------------------------------------------------- 360 #[cfg(test)] 361 mod tests { 362 use super::*; 363 364 #[test] 365 fn example_2_1d() { 366 let options = opt::Opt { 367 ceed_spec: "/cpu/self/ref/serial".to_string(), 368 dim: 1, 369 mesh_degree: 4, 370 solution_degree: 4, 371 num_qpts: 6, 372 problem_size_requested: -1, 373 test: true, 374 quiet: true, 375 gallery: false, 376 }; 377 assert!(example_2(options).is_ok()); 378 } 379 380 #[test] 381 fn example_2_2d() { 382 let options = opt::Opt { 383 ceed_spec: "/cpu/self/ref/serial".to_string(), 384 dim: 2, 385 mesh_degree: 4, 386 solution_degree: 4, 387 num_qpts: 6, 388 problem_size_requested: -1, 389 test: true, 390 quiet: true, 391 gallery: false, 392 }; 393 assert!(example_2(options).is_ok()); 394 } 395 396 #[test] 397 fn example_2_3d() { 398 let options = opt::Opt { 399 ceed_spec: "/cpu/self/ref/serial".to_string(), 400 dim: 3, 401 mesh_degree: 4, 402 solution_degree: 4, 403 num_qpts: 6, 404 problem_size_requested: -1, 405 test: true, 406 quiet: false, 407 gallery: false, 408 }; 409 assert!(example_2(options).is_ok()); 410 } 411 412 #[test] 413 fn example_2_1d_gallery() { 414 let options = opt::Opt { 415 ceed_spec: "/cpu/self/ref/serial".to_string(), 416 dim: 1, 417 mesh_degree: 4, 418 solution_degree: 4, 419 num_qpts: 6, 420 problem_size_requested: -1, 421 test: true, 422 quiet: true, 423 gallery: true, 424 }; 425 assert!(example_2(options).is_ok()); 426 } 427 428 #[test] 429 fn example_2_2d_gallery() { 430 let options = opt::Opt { 431 ceed_spec: "/cpu/self/ref/serial".to_string(), 432 dim: 2, 433 mesh_degree: 4, 434 solution_degree: 4, 435 num_qpts: 6, 436 problem_size_requested: -1, 437 test: true, 438 quiet: true, 439 gallery: true, 440 }; 441 assert!(example_2(options).is_ok()); 442 } 443 444 #[test] 445 fn example_2_3d_gallery() { 446 let options = opt::Opt { 447 ceed_spec: "/cpu/self/ref/serial".to_string(), 448 dim: 3, 449 mesh_degree: 4, 450 solution_degree: 4, 451 num_qpts: 6, 452 problem_size_requested: -1, 453 test: true, 454 quiet: true, 455 gallery: true, 456 }; 457 assert!(example_2(options).is_ok()); 458 } 459 } 460 461 // ---------------------------------------------------------------------------- 462