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() -> libceed::Result<()> { 43 let options = opt::Opt::from_args(); 44 example_2(options) 45 } 46 47 fn example_2(options: opt::Opt) -> libceed::Result<()> { 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 .check()?; 246 247 // Compute the quadrature data for the diff operator 248 let elem_qpts = num_qpts.pow(dim as u32); 249 let num_elem: usize = num_xyz.iter().take(dim).product(); 250 let mut qdata = ceed.vector(num_elem * elem_qpts * dim * (dim + 1) / 2)?; 251 op_build.apply(&mesh_coords, &mut qdata)?; 252 253 // QFunction that applies the diff operator 254 // -- QFunction from user closure 255 let apply_diff = move |[ug, qdata, ..]: QFunctionInputs, [vg, ..]: QFunctionOutputs| { 256 // Apply diffusion operator 257 match dim { 258 1 => vg 259 .iter_mut() 260 .zip(ug.iter().zip(qdata.iter())) 261 .for_each(|(vg, (ug, w))| *vg = ug * w), 262 2 => { 263 let q = qdata.len() / 3; 264 for i in 0..q { 265 let du = [ug[i + q * 0], ug[i + q * 1]]; 266 let dxdxdxdx_t = [ 267 [qdata[i + 0 * q], qdata[i + 2 * q]], 268 [qdata[i + 2 * q], qdata[i + 1 * q]], 269 ]; 270 for j in 0..2 { 271 vg[i + j * q] = du[0] * dxdxdxdx_t[0][j] + du[1] * dxdxdxdx_t[1][j]; 272 } 273 } 274 } 275 3 => { 276 let q = qdata.len() / 6; 277 for i in 0..q { 278 let du = [ug[i + q * 0], ug[i + q * 1], ug[i + q * 2]]; 279 let dxdxdxdx_t = [ 280 [qdata[i + 0 * q], qdata[i + 5 * q], qdata[i + 4 * q]], 281 [qdata[i + 5 * q], qdata[i + 1 * q], qdata[i + 3 * q]], 282 [qdata[i + 4 * q], qdata[i + 3 * q], qdata[i + 2 * q]], 283 ]; 284 for j in 0..3 { 285 vg[i + j * q] = du[0] * dxdxdxdx_t[0][j] 286 + du[1] * dxdxdxdx_t[1][j] 287 + du[2] * dxdxdxdx_t[2][j]; 288 } 289 } 290 } 291 _ => unreachable!(), 292 }; 293 294 // Return clean error code 295 0 296 }; 297 let qf_diff_closure = ceed 298 .q_function_interior(1, Box::new(apply_diff))? 299 .input("du", dim, EvalMode::Grad)? 300 .input("qdata", dim * (dim + 1) / 2, EvalMode::None)? 301 .output("dv", dim, EvalMode::Grad)?; 302 // -- QFunction from gallery 303 let qf_diff_named = { 304 let name = format!("Poisson{}DApply", dim); 305 ceed.q_function_interior_by_name(&name)? 306 }; 307 // -- QFunction for use with Operator 308 let qf_diff = if gallery { 309 QFunctionOpt::SomeQFunctionByName(&qf_diff_named) 310 } else { 311 QFunctionOpt::SomeQFunction(&qf_diff_closure) 312 }; 313 314 // Diff Operator 315 let op_diff = ceed 316 .operator(qf_diff, QFunctionOpt::None, QFunctionOpt::None)? 317 .field("du", &restr_solution, &basis_solution, VectorOpt::Active)? 318 .field("qdata", &restr_qdata, BasisOpt::Collocated, &qdata)? 319 .field("dv", &restr_solution, &basis_solution, VectorOpt::Active)? 320 .check()?; 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: Scalar = 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 => 10000.0 * libceed::EPSILON, 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::Error { 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