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