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