1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative 16 17 //! A Ceed QFunction represents the spatial terms of the point-wise functions 18 //! describing the physics at the quadrature points. 19 20 use std::pin::Pin; 21 22 use crate::prelude::*; 23 24 pub type QFunctionInputs<'a> = [&'a [f64]; MAX_QFUNCTION_FIELDS]; 25 pub type QFunctionOutputs<'a> = [&'a mut [f64]; MAX_QFUNCTION_FIELDS]; 26 27 // ----------------------------------------------------------------------------- 28 // CeedQFunction option 29 // ----------------------------------------------------------------------------- 30 #[derive(Clone, Copy)] 31 pub enum QFunctionOpt<'a> { 32 SomeQFunction(&'a QFunction<'a>), 33 SomeQFunctionByName(&'a QFunctionByName<'a>), 34 None, 35 } 36 37 /// Construct a QFunctionOpt reference from a QFunction reference 38 impl<'a> From<&'a QFunction<'_>> for QFunctionOpt<'a> { 39 fn from(qfunc: &'a QFunction) -> Self { 40 debug_assert!(qfunc.qf_core.ptr != unsafe { bind_ceed::CEED_QFUNCTION_NONE }); 41 Self::SomeQFunction(qfunc) 42 } 43 } 44 45 /// Construct a QFunctionOpt reference from a QFunction by Name reference 46 impl<'a> From<&'a QFunctionByName<'_>> for QFunctionOpt<'a> { 47 fn from(qfunc: &'a QFunctionByName) -> Self { 48 debug_assert!(qfunc.qf_core.ptr != unsafe { bind_ceed::CEED_QFUNCTION_NONE }); 49 Self::SomeQFunctionByName(qfunc) 50 } 51 } 52 53 impl<'a> QFunctionOpt<'a> { 54 /// Transform a Rust libCEED QFunctionOpt into C libCEED CeedQFunction 55 pub(crate) fn to_raw(self) -> bind_ceed::CeedQFunction { 56 match self { 57 Self::SomeQFunction(qfunc) => qfunc.qf_core.ptr, 58 Self::SomeQFunctionByName(qfunc) => qfunc.qf_core.ptr, 59 Self::None => unsafe { bind_ceed::CEED_QFUNCTION_NONE }, 60 } 61 } 62 } 63 64 // ----------------------------------------------------------------------------- 65 // CeedQFunction context wrapper 66 // ----------------------------------------------------------------------------- 67 pub(crate) struct QFunctionCore<'a> { 68 ceed: &'a crate::Ceed, 69 ptr: bind_ceed::CeedQFunction, 70 } 71 72 struct QFunctionTrampolineData { 73 number_inputs: usize, 74 number_outputs: usize, 75 input_sizes: [usize; MAX_QFUNCTION_FIELDS], 76 output_sizes: [usize; MAX_QFUNCTION_FIELDS], 77 user_f: Box<QFunctionUserClosure>, 78 } 79 80 pub struct QFunction<'a> { 81 qf_core: QFunctionCore<'a>, 82 qf_ctx_ptr: bind_ceed::CeedQFunctionContext, 83 trampoline_data: Pin<Box<QFunctionTrampolineData>>, 84 } 85 86 pub struct QFunctionByName<'a> { 87 qf_core: QFunctionCore<'a>, 88 } 89 90 // ----------------------------------------------------------------------------- 91 // Destructor 92 // ----------------------------------------------------------------------------- 93 impl<'a> Drop for QFunctionCore<'a> { 94 fn drop(&mut self) { 95 unsafe { 96 if self.ptr != bind_ceed::CEED_QFUNCTION_NONE { 97 bind_ceed::CeedQFunctionDestroy(&mut self.ptr); 98 } 99 } 100 } 101 } 102 103 impl<'a> Drop for QFunction<'a> { 104 fn drop(&mut self) { 105 unsafe { 106 bind_ceed::CeedQFunctionContextDestroy(&mut self.qf_ctx_ptr); 107 } 108 } 109 } 110 111 // ----------------------------------------------------------------------------- 112 // Display 113 // ----------------------------------------------------------------------------- 114 impl<'a> fmt::Display for QFunctionCore<'a> { 115 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { 116 let mut ptr = std::ptr::null_mut(); 117 let mut sizeloc = crate::MAX_BUFFER_LENGTH; 118 let cstring = unsafe { 119 let file = bind_ceed::open_memstream(&mut ptr, &mut sizeloc); 120 bind_ceed::CeedQFunctionView(self.ptr, file); 121 bind_ceed::fclose(file); 122 CString::from_raw(ptr) 123 }; 124 cstring.to_string_lossy().fmt(f) 125 } 126 } 127 /// View a QFunction 128 /// 129 /// ``` 130 /// # use libceed::prelude::*; 131 /// # let ceed = libceed::Ceed::default_init(); 132 /// let mut user_f = |[u, weights, ..]: QFunctionInputs, [v, ..]: QFunctionOutputs| { 133 /// // Iterate over quadrature points 134 /// v.iter_mut() 135 /// .zip(u.iter().zip(weights.iter())) 136 /// .for_each(|(v, (u, w))| *v = u * w); 137 /// 138 /// // Return clean error code 139 /// 0 140 /// }; 141 /// 142 /// let qf = ceed 143 /// .q_function_interior(1, Box::new(user_f)) 144 /// .unwrap() 145 /// .input("u", 1, EvalMode::Interp) 146 /// .unwrap() 147 /// .input("weights", 1, EvalMode::Weight) 148 /// .unwrap() 149 /// .output("v", 1, EvalMode::Interp) 150 /// .unwrap(); 151 /// 152 /// println!("{}", qf); 153 /// ``` 154 impl<'a> fmt::Display for QFunction<'a> { 155 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { 156 self.qf_core.fmt(f) 157 } 158 } 159 160 /// View a QFunction by Name 161 /// 162 /// ``` 163 /// # use libceed::prelude::*; 164 /// # let ceed = libceed::Ceed::default_init(); 165 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild").unwrap(); 166 /// println!("{}", qf); 167 /// ``` 168 impl<'a> fmt::Display for QFunctionByName<'a> { 169 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { 170 self.qf_core.fmt(f) 171 } 172 } 173 174 // ----------------------------------------------------------------------------- 175 // Core functionality 176 // ----------------------------------------------------------------------------- 177 impl<'a> QFunctionCore<'a> { 178 // Common implementation 179 pub fn apply(&self, Q: usize, u: &[Vector], v: &[Vector]) -> crate::Result<i32> { 180 let mut u_c = [std::ptr::null_mut(); MAX_QFUNCTION_FIELDS]; 181 for i in 0..std::cmp::min(MAX_QFUNCTION_FIELDS, u.len()) { 182 u_c[i] = u[i].ptr; 183 } 184 let mut v_c = [std::ptr::null_mut(); MAX_QFUNCTION_FIELDS]; 185 for i in 0..std::cmp::min(MAX_QFUNCTION_FIELDS, v.len()) { 186 v_c[i] = v[i].ptr; 187 } 188 let Q = i32::try_from(Q).unwrap(); 189 let ierr = unsafe { 190 bind_ceed::CeedQFunctionApply(self.ptr, Q, u_c.as_mut_ptr(), v_c.as_mut_ptr()) 191 }; 192 self.ceed.check_error(ierr) 193 } 194 } 195 196 // ----------------------------------------------------------------------------- 197 // User QFunction Closure 198 // ----------------------------------------------------------------------------- 199 pub type QFunctionUserClosure = 200 dyn FnMut([&[f64]; MAX_QFUNCTION_FIELDS], [&mut [f64]; MAX_QFUNCTION_FIELDS]) -> i32; 201 202 macro_rules! mut_max_fields { 203 ($e:expr) => { 204 [ 205 $e, $e, $e, $e, $e, $e, $e, $e, $e, $e, $e, $e, $e, $e, $e, $e, 206 ] 207 }; 208 } 209 unsafe extern "C" fn trampoline( 210 ctx: *mut ::std::os::raw::c_void, 211 q: bind_ceed::CeedInt, 212 inputs: *const *const bind_ceed::CeedScalar, 213 outputs: *const *mut bind_ceed::CeedScalar, 214 ) -> ::std::os::raw::c_int { 215 let trampoline_data: Pin<&mut QFunctionTrampolineData> = std::mem::transmute(ctx); 216 217 // Inputs 218 let inputs_slice: &[*const bind_ceed::CeedScalar] = 219 std::slice::from_raw_parts(inputs, MAX_QFUNCTION_FIELDS); 220 let mut inputs_array: [&[f64]; MAX_QFUNCTION_FIELDS] = [&[0.0]; MAX_QFUNCTION_FIELDS]; 221 inputs_slice 222 .iter() 223 .enumerate() 224 .map(|(i, &x)| { 225 std::slice::from_raw_parts(x, trampoline_data.input_sizes[i] * q as usize) as &[f64] 226 }) 227 .zip(inputs_array.iter_mut()) 228 .for_each(|(x, a)| *a = x); 229 230 // Outputs 231 let outputs_slice: &[*mut bind_ceed::CeedScalar] = 232 std::slice::from_raw_parts(outputs, MAX_QFUNCTION_FIELDS); 233 let mut outputs_array: [&mut [f64]; MAX_QFUNCTION_FIELDS] = mut_max_fields!(&mut [0.0]); 234 outputs_slice 235 .iter() 236 .enumerate() 237 .map(|(i, &x)| { 238 std::slice::from_raw_parts_mut(x, trampoline_data.output_sizes[i] * q as usize) 239 as &mut [f64] 240 }) 241 .zip(outputs_array.iter_mut()) 242 .for_each(|(x, a)| *a = x); 243 244 // User closure 245 (trampoline_data.get_unchecked_mut().user_f)(inputs_array, outputs_array) 246 } 247 248 // ----------------------------------------------------------------------------- 249 // QFunction 250 // ----------------------------------------------------------------------------- 251 impl<'a> QFunction<'a> { 252 // Constructor 253 pub fn create( 254 ceed: &'a crate::Ceed, 255 vlength: usize, 256 user_f: Box<QFunctionUserClosure>, 257 ) -> crate::Result<Self> { 258 let source_c = CString::new("").expect("CString::new failed"); 259 let mut ptr = std::ptr::null_mut(); 260 261 // Context for closure 262 let number_inputs = 0; 263 let number_outputs = 0; 264 let input_sizes = [0; MAX_QFUNCTION_FIELDS]; 265 let output_sizes = [0; MAX_QFUNCTION_FIELDS]; 266 let trampoline_data = unsafe { 267 Pin::new_unchecked(Box::new(QFunctionTrampolineData { 268 number_inputs, 269 number_outputs, 270 input_sizes, 271 output_sizes, 272 user_f, 273 })) 274 }; 275 276 // Create QFunction 277 let vlength = i32::try_from(vlength).unwrap(); 278 let mut ierr = unsafe { 279 bind_ceed::CeedQFunctionCreateInterior( 280 ceed.ptr, 281 vlength, 282 Some(trampoline), 283 source_c.as_ptr(), 284 &mut ptr, 285 ) 286 }; 287 ceed.check_error(ierr)?; 288 289 // Set closure 290 let mut qf_ctx_ptr = std::ptr::null_mut(); 291 ierr = unsafe { bind_ceed::CeedQFunctionContextCreate(ceed.ptr, &mut qf_ctx_ptr) }; 292 ceed.check_error(ierr)?; 293 ierr = unsafe { 294 bind_ceed::CeedQFunctionContextSetData( 295 qf_ctx_ptr, 296 crate::MemType::Host as bind_ceed::CeedMemType, 297 crate::CopyMode::UsePointer as bind_ceed::CeedCopyMode, 298 std::mem::size_of::<QFunctionTrampolineData>() as u64, 299 std::mem::transmute(trampoline_data.as_ref()), 300 ) 301 }; 302 ceed.check_error(ierr)?; 303 ierr = unsafe { bind_ceed::CeedQFunctionSetContext(ptr, qf_ctx_ptr) }; 304 ceed.check_error(ierr)?; 305 Ok(Self { 306 qf_core: QFunctionCore { ceed, ptr }, 307 qf_ctx_ptr, 308 trampoline_data, 309 }) 310 } 311 312 /// Apply the action of a QFunction 313 /// 314 /// * `Q` - The number of quadrature points 315 /// * `input` - Array of input Vectors 316 /// * `output` - Array of output Vectors 317 /// 318 /// ``` 319 /// # use libceed::prelude::*; 320 /// # let ceed = libceed::Ceed::default_init(); 321 /// let mut user_f = |[u, weights, ..]: QFunctionInputs, [v, ..]: QFunctionOutputs| { 322 /// // Iterate over quadrature points 323 /// v.iter_mut() 324 /// .zip(u.iter().zip(weights.iter())) 325 /// .for_each(|(v, (u, w))| *v = u * w); 326 /// 327 /// // Return clean error code 328 /// 0 329 /// }; 330 /// 331 /// let qf = ceed 332 /// .q_function_interior(1, Box::new(user_f)) 333 /// .unwrap() 334 /// .input("u", 1, EvalMode::Interp) 335 /// .unwrap() 336 /// .input("weights", 1, EvalMode::Weight) 337 /// .unwrap() 338 /// .output("v", 1, EvalMode::Interp) 339 /// .unwrap(); 340 /// 341 /// const Q: usize = 8; 342 /// let mut w = [0.; Q]; 343 /// let mut u = [0.; Q]; 344 /// let mut v = [0.; Q]; 345 /// 346 /// for i in 0..Q { 347 /// let x = 2. * (i as f64) / ((Q as f64) - 1.) - 1.; 348 /// u[i] = 2. + 3. * x + 5. * x * x; 349 /// w[i] = 1. - x * x; 350 /// v[i] = u[i] * w[i]; 351 /// } 352 /// 353 /// let uu = ceed.vector_from_slice(&u).unwrap(); 354 /// let ww = ceed.vector_from_slice(&w).unwrap(); 355 /// let mut vv = ceed.vector(Q).unwrap(); 356 /// vv.set_value(0.0); 357 /// { 358 /// let input = vec![uu, ww]; 359 /// let mut output = vec![vv]; 360 /// qf.apply(Q, &input, &output).unwrap(); 361 /// vv = output.remove(0); 362 /// } 363 /// 364 /// vv.view() 365 /// .iter() 366 /// .zip(v.iter()) 367 /// .for_each(|(computed, actual)| { 368 /// assert_eq!( 369 /// *computed, *actual, 370 /// "Incorrect value in QFunction application" 371 /// ); 372 /// }); 373 /// ``` 374 pub fn apply(&self, Q: usize, u: &[Vector], v: &[Vector]) -> crate::Result<i32> { 375 self.qf_core.apply(Q, u, v) 376 } 377 378 /// Add a QFunction input 379 /// 380 /// * `fieldname` - Name of QFunction field 381 /// * `size` - Size of QFunction field, `(ncomp * dim)` for `Grad` or 382 /// `(ncomp * 1)` for `None`, `Interp`, and `Weight` 383 /// * `emode` - `EvalMode::None` to use values directly, `EvalMode::Interp` 384 /// to use interpolated values, `EvalMode::Grad` to use 385 /// gradients, `EvalMode::Weight` to use quadrature weights 386 /// 387 /// ``` 388 /// # use libceed::prelude::*; 389 /// # let ceed = libceed::Ceed::default_init(); 390 /// let mut user_f = |[u, weights, ..]: QFunctionInputs, [v, ..]: QFunctionOutputs| { 391 /// // Iterate over quadrature points 392 /// v.iter_mut() 393 /// .zip(u.iter().zip(weights.iter())) 394 /// .for_each(|(v, (u, w))| *v = u * w); 395 /// 396 /// // Return clean error code 397 /// 0 398 /// }; 399 /// 400 /// let mut qf = ceed.q_function_interior(1, Box::new(user_f)).unwrap(); 401 /// 402 /// qf = qf.input("u", 1, EvalMode::Interp).unwrap(); 403 /// qf = qf.input("weights", 1, EvalMode::Weight).unwrap(); 404 /// ``` 405 pub fn input( 406 mut self, 407 fieldname: &str, 408 size: usize, 409 emode: crate::EvalMode, 410 ) -> crate::Result<Self> { 411 let name_c = CString::new(fieldname).expect("CString::new failed"); 412 let idx = self.trampoline_data.number_inputs; 413 self.trampoline_data.input_sizes[idx] = size; 414 self.trampoline_data.number_inputs += 1; 415 let (size, emode) = ( 416 i32::try_from(size).unwrap(), 417 emode as bind_ceed::CeedEvalMode, 418 ); 419 let ierr = unsafe { 420 bind_ceed::CeedQFunctionAddInput(self.qf_core.ptr, name_c.as_ptr(), size, emode) 421 }; 422 self.qf_core.ceed.check_error(ierr)?; 423 Ok(self) 424 } 425 426 /// Add a QFunction output 427 /// 428 /// * `fieldname` - Name of QFunction field 429 /// * `size` - Size of QFunction field, `(ncomp * dim)` for `Grad` or 430 /// `(ncomp * 1)` for `None` and `Interp` 431 /// * `emode` - `EvalMode::None` to use values directly, `EvalMode::Interp` 432 /// to use interpolated values, `EvalMode::Grad` to use 433 /// gradients 434 /// 435 /// ``` 436 /// # use libceed::prelude::*; 437 /// # let ceed = libceed::Ceed::default_init(); 438 /// let mut user_f = |[u, weights, ..]: QFunctionInputs, [v, ..]: QFunctionOutputs| { 439 /// // Iterate over quadrature points 440 /// v.iter_mut() 441 /// .zip(u.iter().zip(weights.iter())) 442 /// .for_each(|(v, (u, w))| *v = u * w); 443 /// 444 /// // Return clean error code 445 /// 0 446 /// }; 447 /// 448 /// let mut qf = ceed.q_function_interior(1, Box::new(user_f)).unwrap(); 449 /// 450 /// qf.output("v", 1, EvalMode::Interp).unwrap(); 451 /// ``` 452 pub fn output( 453 mut self, 454 fieldname: &str, 455 size: usize, 456 emode: crate::EvalMode, 457 ) -> crate::Result<Self> { 458 let name_c = CString::new(fieldname).expect("CString::new failed"); 459 let idx = self.trampoline_data.number_outputs; 460 self.trampoline_data.output_sizes[idx] = size; 461 self.trampoline_data.number_outputs += 1; 462 let (size, emode) = ( 463 i32::try_from(size).unwrap(), 464 emode as bind_ceed::CeedEvalMode, 465 ); 466 let ierr = unsafe { 467 bind_ceed::CeedQFunctionAddOutput(self.qf_core.ptr, name_c.as_ptr(), size, emode) 468 }; 469 self.qf_core.ceed.check_error(ierr)?; 470 Ok(self) 471 } 472 } 473 474 // ----------------------------------------------------------------------------- 475 // QFunction 476 // ----------------------------------------------------------------------------- 477 impl<'a> QFunctionByName<'a> { 478 // Constructor 479 pub fn create(ceed: &'a crate::Ceed, name: &str) -> crate::Result<Self> { 480 let name_c = CString::new(name).expect("CString::new failed"); 481 let mut ptr = std::ptr::null_mut(); 482 let ierr = unsafe { 483 bind_ceed::CeedQFunctionCreateInteriorByName(ceed.ptr, name_c.as_ptr(), &mut ptr) 484 }; 485 ceed.check_error(ierr)?; 486 Ok(Self { 487 qf_core: QFunctionCore { ceed, ptr }, 488 }) 489 } 490 491 /// Apply the action of a QFunction 492 /// 493 /// * `Q` - The number of quadrature points 494 /// * `input` - Array of input Vectors 495 /// * `output` - Array of output Vectors 496 /// 497 /// ``` 498 /// # use libceed::prelude::*; 499 /// # let ceed = libceed::Ceed::default_init(); 500 /// const Q: usize = 8; 501 /// let qf_build = ceed.q_function_interior_by_name("Mass1DBuild").unwrap(); 502 /// let qf_mass = ceed.q_function_interior_by_name("MassApply").unwrap(); 503 /// 504 /// let mut j = [0.; Q]; 505 /// let mut w = [0.; Q]; 506 /// let mut u = [0.; Q]; 507 /// let mut v = [0.; Q]; 508 /// 509 /// for i in 0..Q { 510 /// let x = 2. * (i as f64) / ((Q as f64) - 1.) - 1.; 511 /// j[i] = 1.; 512 /// w[i] = 1. - x * x; 513 /// u[i] = 2. + 3. * x + 5. * x * x; 514 /// v[i] = w[i] * u[i]; 515 /// } 516 /// 517 /// let jj = ceed.vector_from_slice(&j).unwrap(); 518 /// let ww = ceed.vector_from_slice(&w).unwrap(); 519 /// let uu = ceed.vector_from_slice(&u).unwrap(); 520 /// let mut vv = ceed.vector(Q).unwrap(); 521 /// vv.set_value(0.0); 522 /// let mut qdata = ceed.vector(Q).unwrap(); 523 /// qdata.set_value(0.0); 524 /// 525 /// { 526 /// let mut input = vec![jj, ww]; 527 /// let mut output = vec![qdata]; 528 /// qf_build.apply(Q, &input, &output).unwrap(); 529 /// qdata = output.remove(0); 530 /// } 531 /// 532 /// { 533 /// let mut input = vec![qdata, uu]; 534 /// let mut output = vec![vv]; 535 /// qf_mass.apply(Q, &input, &output).unwrap(); 536 /// vv = output.remove(0); 537 /// } 538 /// 539 /// vv.view() 540 /// .iter() 541 /// .zip(v.iter()) 542 /// .for_each(|(computed, actual)| { 543 /// assert_eq!( 544 /// *computed, *actual, 545 /// "Incorrect value in QFunction application" 546 /// ); 547 /// }); 548 /// ``` 549 pub fn apply(&self, Q: usize, u: &[Vector], v: &[Vector]) -> crate::Result<i32> { 550 self.qf_core.apply(Q, u, v) 551 } 552 } 553 554 // ----------------------------------------------------------------------------- 555