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 (restr_mesh, _) = 107 mesh::build_cartesian_restriction(&ceed, dim, num_xyz, mesh_degree, ncomp_x, num_qpts)?; 108 let (_, restr_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 (restr_solution, _) = 118 mesh::build_cartesian_restriction(&ceed, dim, num_xyz, solution_degree, 1, num_qpts)?; 119 let mesh_size = restr_mesh.lvector_size(); 120 let solution_size = restr_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", &restr_mesh, &basis_mesh, VectorOpt::Active)? 225 .field( 226 "weights", 227 ElemRestrictionOpt::None, 228 &basis_mesh, 229 VectorOpt::None, 230 )? 231 .field( 232 "qdata", 233 &restr_qdata, 234 BasisOpt::Collocated, 235 VectorOpt::Active, 236 )? 237 .check()?; 238 239 // Compute the quadrature data for the diff operator 240 let elem_qpts = num_qpts.pow(dim as u32); 241 let num_elem: usize = num_xyz.iter().take(dim).product(); 242 let mut qdata = ceed.vector(num_elem * elem_qpts * dim * (dim + 1) / 2)?; 243 op_build.apply(&mesh_coords, &mut qdata)?; 244 245 // QFunction that applies the diff operator 246 // -- QFunction from user closure 247 let apply_diff = move |[ug, qdata, ..]: QFunctionInputs, [vg, ..]: QFunctionOutputs| { 248 // Apply diffusion operator 249 match dim { 250 1 => vg 251 .iter_mut() 252 .zip(ug.iter().zip(qdata.iter())) 253 .for_each(|(vg, (ug, w))| *vg = ug * w), 254 2 => { 255 let q = qdata.len() / 3; 256 for i in 0..q { 257 let du = [ug[i + q * 0], ug[i + q * 1]]; 258 let dxdxdxdx_t = [ 259 [qdata[i + 0 * q], qdata[i + 2 * q]], 260 [qdata[i + 2 * q], qdata[i + 1 * q]], 261 ]; 262 for j in 0..2 { 263 vg[i + j * q] = du[0] * dxdxdxdx_t[0][j] + du[1] * dxdxdxdx_t[1][j]; 264 } 265 } 266 } 267 3 => { 268 let q = qdata.len() / 6; 269 for i in 0..q { 270 let du = [ug[i + q * 0], ug[i + q * 1], ug[i + q * 2]]; 271 let dxdxdxdx_t = [ 272 [qdata[i + 0 * q], qdata[i + 5 * q], qdata[i + 4 * q]], 273 [qdata[i + 5 * q], qdata[i + 1 * q], qdata[i + 3 * q]], 274 [qdata[i + 4 * q], qdata[i + 3 * q], qdata[i + 2 * q]], 275 ]; 276 for j in 0..3 { 277 vg[i + j * q] = du[0] * dxdxdxdx_t[0][j] 278 + du[1] * dxdxdxdx_t[1][j] 279 + du[2] * dxdxdxdx_t[2][j]; 280 } 281 } 282 } 283 _ => unreachable!(), 284 }; 285 286 // Return clean error code 287 0 288 }; 289 let qf_diff_closure = ceed 290 .q_function_interior(1, Box::new(apply_diff))? 291 .input("du", dim, EvalMode::Grad)? 292 .input("qdata", dim * (dim + 1) / 2, EvalMode::None)? 293 .output("dv", dim, EvalMode::Grad)?; 294 // -- QFunction from gallery 295 let qf_diff_named = { 296 let name = format!("Poisson{}DApply", dim); 297 ceed.q_function_interior_by_name(&name)? 298 }; 299 // -- QFunction for use with Operator 300 let qf_diff = if gallery { 301 QFunctionOpt::SomeQFunctionByName(&qf_diff_named) 302 } else { 303 QFunctionOpt::SomeQFunction(&qf_diff_closure) 304 }; 305 306 // Diff Operator 307 let op_diff = ceed 308 .operator(qf_diff, QFunctionOpt::None, QFunctionOpt::None)? 309 .name("Poisson")? 310 .field("du", &restr_solution, &basis_solution, VectorOpt::Active)? 311 .field("qdata", &restr_qdata, BasisOpt::Collocated, &qdata)? 312 .field("dv", &restr_solution, &basis_solution, VectorOpt::Active)? 313 .check()?; 314 315 // Solution vectors 316 let mut u = ceed.vector(solution_size)?; 317 let mut v = ceed.vector(solution_size)?; 318 319 // Initialize u with sum of node coordinates 320 let coords = mesh_coords.view()?; 321 u.set_value(0.0)?; 322 for (i, u) in u.view_mut()?.iter_mut().enumerate() { 323 *u = (0..dim).map(|d| coords[i + d * solution_size]).sum(); 324 } 325 326 // Apply the diff operator 327 op_diff.apply(&u, &mut v)?; 328 329 // Compute the mesh surface area 330 let area: Scalar = v.view()?.iter().map(|v| (*v).abs()).sum(); 331 332 // Output results 333 if !quiet { 334 println!("Exact mesh surface area : {:.12}", exact_area); 335 println!("Computed mesh surface_area : {:.12}", area); 336 println!("Surface area error : {:.12e}", area - exact_area); 337 } 338 let tolerance = match dim { 339 1 => 10000.0 * libceed::EPSILON, 340 _ => 1E-1, 341 }; 342 let error = (area - exact_area).abs(); 343 if error > tolerance { 344 println!("Volume error too large: {:.12e}", error); 345 return Err(libceed::Error { 346 message: format!( 347 "Volume error too large - expected: {:.12e}, actual: {:.12e}", 348 tolerance, error 349 ), 350 }); 351 } 352 Ok(()) 353 } 354 355 // ---------------------------------------------------------------------------- 356 // Tests 357 // ---------------------------------------------------------------------------- 358 #[cfg(test)] 359 mod tests { 360 use super::*; 361 362 #[test] 363 fn example_2_1d() { 364 let options = opt::Opt { 365 ceed_spec: "/cpu/self/ref/serial".to_string(), 366 dim: 1, 367 mesh_degree: 4, 368 solution_degree: 4, 369 num_qpts: 6, 370 problem_size_requested: -1, 371 test: true, 372 quiet: true, 373 gallery: false, 374 }; 375 assert!(example_2(options).is_ok()); 376 } 377 378 #[test] 379 fn example_2_2d() { 380 let options = opt::Opt { 381 ceed_spec: "/cpu/self/ref/serial".to_string(), 382 dim: 2, 383 mesh_degree: 4, 384 solution_degree: 4, 385 num_qpts: 6, 386 problem_size_requested: -1, 387 test: true, 388 quiet: true, 389 gallery: false, 390 }; 391 assert!(example_2(options).is_ok()); 392 } 393 394 #[test] 395 fn example_2_3d() { 396 let options = opt::Opt { 397 ceed_spec: "/cpu/self/ref/serial".to_string(), 398 dim: 3, 399 mesh_degree: 4, 400 solution_degree: 4, 401 num_qpts: 6, 402 problem_size_requested: -1, 403 test: true, 404 quiet: false, 405 gallery: false, 406 }; 407 assert!(example_2(options).is_ok()); 408 } 409 410 #[test] 411 fn example_2_1d_gallery() { 412 let options = opt::Opt { 413 ceed_spec: "/cpu/self/ref/serial".to_string(), 414 dim: 1, 415 mesh_degree: 4, 416 solution_degree: 4, 417 num_qpts: 6, 418 problem_size_requested: -1, 419 test: true, 420 quiet: true, 421 gallery: true, 422 }; 423 assert!(example_2(options).is_ok()); 424 } 425 426 #[test] 427 fn example_2_2d_gallery() { 428 let options = opt::Opt { 429 ceed_spec: "/cpu/self/ref/serial".to_string(), 430 dim: 2, 431 mesh_degree: 4, 432 solution_degree: 4, 433 num_qpts: 6, 434 problem_size_requested: -1, 435 test: true, 436 quiet: true, 437 gallery: true, 438 }; 439 assert!(example_2(options).is_ok()); 440 } 441 442 #[test] 443 fn example_2_3d_gallery() { 444 let options = opt::Opt { 445 ceed_spec: "/cpu/self/ref/serial".to_string(), 446 dim: 3, 447 mesh_degree: 4, 448 solution_degree: 4, 449 num_qpts: 6, 450 problem_size_requested: -1, 451 test: true, 452 quiet: true, 453 gallery: true, 454 }; 455 assert!(example_2(options).is_ok()); 456 } 457 } 458 459 // ---------------------------------------------------------------------------- 460