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.set_value(0.0)?; 329 for (i, u) in u.view_mut()?.iter_mut().enumerate() { 330 *u = (0..dim).map(|d| coords[i + d * solution_size]).sum(); 331 } 332 333 // Apply the diff operator 334 op_diff.apply(&u, &mut v)?; 335 336 // Compute the mesh surface area 337 let area: Scalar = v.view()?.iter().map(|v| (*v).abs()).sum(); 338 339 // Output results 340 if !quiet { 341 println!("Exact mesh surface area : {:.12}", exact_area); 342 println!("Computed mesh surface_area : {:.12}", area); 343 println!("Surface area error : {:.12e}", area - exact_area); 344 } 345 let tolerance = match dim { 346 1 => 10000.0 * libceed::EPSILON, 347 _ => 1E-1, 348 }; 349 let error = (area - exact_area).abs(); 350 if error > tolerance { 351 println!("Volume error too large: {:.12e}", error); 352 return Err(libceed::Error { 353 message: format!( 354 "Volume error too large - expected: {:.12e}, actual: {:.12e}", 355 tolerance, error 356 ), 357 }); 358 } 359 Ok(()) 360 } 361 362 // ---------------------------------------------------------------------------- 363 // Tests 364 // ---------------------------------------------------------------------------- 365 #[cfg(test)] 366 mod tests { 367 use super::*; 368 369 #[test] 370 fn example_2_1d() { 371 let options = opt::Opt { 372 ceed_spec: "/cpu/self/ref/serial".to_string(), 373 dim: 1, 374 mesh_degree: 4, 375 solution_degree: 4, 376 num_qpts: 6, 377 problem_size_requested: -1, 378 test: true, 379 quiet: true, 380 gallery: false, 381 }; 382 assert!(example_2(options).is_ok()); 383 } 384 385 #[test] 386 fn example_2_2d() { 387 let options = opt::Opt { 388 ceed_spec: "/cpu/self/ref/serial".to_string(), 389 dim: 2, 390 mesh_degree: 4, 391 solution_degree: 4, 392 num_qpts: 6, 393 problem_size_requested: -1, 394 test: true, 395 quiet: true, 396 gallery: false, 397 }; 398 assert!(example_2(options).is_ok()); 399 } 400 401 #[test] 402 fn example_2_3d() { 403 let options = opt::Opt { 404 ceed_spec: "/cpu/self/ref/serial".to_string(), 405 dim: 3, 406 mesh_degree: 4, 407 solution_degree: 4, 408 num_qpts: 6, 409 problem_size_requested: -1, 410 test: true, 411 quiet: false, 412 gallery: false, 413 }; 414 assert!(example_2(options).is_ok()); 415 } 416 417 #[test] 418 fn example_2_1d_gallery() { 419 let options = opt::Opt { 420 ceed_spec: "/cpu/self/ref/serial".to_string(), 421 dim: 1, 422 mesh_degree: 4, 423 solution_degree: 4, 424 num_qpts: 6, 425 problem_size_requested: -1, 426 test: true, 427 quiet: true, 428 gallery: true, 429 }; 430 assert!(example_2(options).is_ok()); 431 } 432 433 #[test] 434 fn example_2_2d_gallery() { 435 let options = opt::Opt { 436 ceed_spec: "/cpu/self/ref/serial".to_string(), 437 dim: 2, 438 mesh_degree: 4, 439 solution_degree: 4, 440 num_qpts: 6, 441 problem_size_requested: -1, 442 test: true, 443 quiet: true, 444 gallery: true, 445 }; 446 assert!(example_2(options).is_ok()); 447 } 448 449 #[test] 450 fn example_2_3d_gallery() { 451 let options = opt::Opt { 452 ceed_spec: "/cpu/self/ref/serial".to_string(), 453 dim: 3, 454 mesh_degree: 4, 455 solution_degree: 4, 456 num_qpts: 6, 457 problem_size_requested: -1, 458 test: true, 459 quiet: true, 460 gallery: true, 461 }; 462 assert!(example_2(options).is_ok()); 463 } 464 } 465 466 // ---------------------------------------------------------------------------- 467