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