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