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