1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254
//! A wrapper around cvode and cvodes from the sundials tool suite. //! //! Users should be mostly interested in [`SolverSensi`] and [`SolverNoSensi`]. //! //! # Building sundials //! //! To build sundials, activate the `sundials-sys/build_libraries` feature. //! //! # Examples //! //! ## Oscillator //! //! An oscillatory system defined by `x'' = -k * x`. //! //! ### Without sensitivities //! //! ```rust //! use cvode_wrap::*; //! let y0 = [0., 1.]; //! //define the right-hand-side //! fn f(_t: Realtype, y: &[Realtype; 2], ydot: &mut [Realtype; 2], k: &Realtype) -> RhsResult { //! *ydot = [y[1], -y[0] * k]; //! RhsResult::Ok //! } //! //initialize the solver //! let mut solver = SolverNoSensi::new( //! LinearMultistepMethod::Adams, //! f, //! 0., //! &y0, //! 1e-4, //! AbsTolerance::scalar(1e-4), //! 1e-2, //! ) //! .unwrap(); //! //and solve //! let ts: Vec<_> = (1..100).collect(); //! println!("0,{},{}", y0[0], y0[1]); //! for &t in &ts { //! let (_tret, &[x, xdot]) = solver.step(t as _, StepKind::Normal).unwrap(); //! println!("{},{},{}", t, x, xdot); //! } //! ``` //! //! ### With sensitivities //! //! The sensitivities are computed with respect to `x(0)`, `x'(0)` and `k`. //! //! ```rust //! use cvode_wrap::*; //! let y0 = [0., 1.]; //! //define the right-hand-side //! fn f(_t: Realtype, y: &[Realtype; 2], ydot: &mut [Realtype; 2], k: &Realtype) -> RhsResult { //! *ydot = [y[1], -y[0] * k]; //! RhsResult::Ok //! } //! //define the sensitivity function for the right hand side //! fn fs( //! _t: Realtype, //! y: &[Realtype; 2], //! _ydot: &[Realtype; 2], //! ys: [&[Realtype; 2]; N_SENSI], //! ysdot: [&mut [Realtype; 2]; N_SENSI], //! k: &Realtype, //! ) -> RhsResult { //! // Mind that when indexing sensitivities, the first index //! // is the parameter index, and the second the state variable //! // index //! *ysdot[0] = [ys[0][1], -ys[0][0] * k]; //! *ysdot[1] = [ys[1][1], -ys[1][0] * k]; //! *ysdot[2] = [ys[2][1], -ys[2][0] * k - y[0]]; //! RhsResult::Ok //! } //! //! const N_SENSI: usize = 3; //! //! // the sensitivities in order are d/dy0[0], d/dy0[1] and d/dk //! let ys0 = [[1., 0.], [0., 1.], [0., 0.]]; //! //! //initialize the solver //! let mut solver = SolverSensi::new( //! LinearMultistepMethod::Adams, //! f, //! fs, //! 0., //! &y0, //! &ys0, //! 1e-4, //! AbsTolerance::scalar(1e-4), //! SensiAbsTolerance::scalar([1e-4; N_SENSI]), //! 1e-2, //! ) //! .unwrap(); //! //and solve //! let ts: Vec<_> = (1..100).collect(); //! println!("0,{},{}", y0[0], y0[1]); //! for &t in &ts { //! let (_tret, &[x, xdot], [&[dy0_dy00, dy1_dy00], &[dy0_dy01, dy1_dy01], &[dy0_dk, dy1_dk]]) = //! solver.step(t as _, StepKind::Normal).unwrap(); //! println!( //! "{},{},{},{},{},{},{},{},{}", //! t, x, xdot, dy0_dy00, dy1_dy00, dy0_dy01, dy1_dy01, dy0_dk, dy1_dk //! ); //! } //! ``` use std::{ffi::c_void, os::raw::c_int, ptr::NonNull}; use sundials_sys::realtype; mod nvector; pub use nvector::{NVectorSerial, NVectorSerialHeapAllocated}; mod cvode; mod cvode_sens; pub use cvode::Solver as SolverNoSensi; pub use cvode_sens::Solver as SolverSensi; /// The floatting-point type sundials was compiled with pub type Realtype = realtype; #[repr(i32)] #[derive(Debug)] /// An integration method. pub enum LinearMultistepMethod { /// Recomended for non-stiff problems. Adams = sundials_sys::CV_ADAMS, /// Recommended for stiff problems. Bdf = sundials_sys::CV_BDF, } /// A return type for the right-hand-side rust function. /// /// Adapted from Sundials cv-ode guide version 5.7 (BSD Licensed), setcion 4.6.1 : /// /// > If a recoverable error occurred, `cvode` will attempt to correct, /// > if the error is unrecoverable, the integration is halted. /// > /// > A recoverable failure error return is typically used to flag a value of /// > the dependent variableythat is “illegal” in some way (e.g., negative where /// > only a non-negative value is physically meaningful). If such a return is /// > made, `cvode` will attempt to recover (possibly repeating the nonlinear solve, /// > or reducing the step size) in order to avoid this recoverable error return. pub enum RhsResult { /// Indicates that there was no error Ok, /// Indicate that there was a recoverable error and its code RecoverableError(u8), /// Indicatest hat there was a non recoverable error NonRecoverableError(u8), } /// Type of integration step #[repr(i32)] pub enum StepKind { /// The `NORMAL`option causes the solver to take internal steps /// until it has reached or just passed the user-specified time. /// The solver then interpolates in order to return an approximate /// value of y at the desired time. Normal = sundials_sys::CV_NORMAL, /// The `CV_ONE_STEP` option tells the solver to take just one /// internal step and then return thesolution at the point reached /// by that step. OneStep = sundials_sys::CV_ONE_STEP, } /// The error type for this crate #[derive(Debug)] pub enum Error { NullPointerError { func_id: &'static str }, ErrorCode { func_id: &'static str, flag: c_int }, } /// An enum representing the choice between a scalar or vector absolute tolerance pub enum AbsTolerance<const SIZE: usize> { Scalar(Realtype), Vector(NVectorSerialHeapAllocated<SIZE>), } impl<const SIZE: usize> AbsTolerance<SIZE> { pub fn scalar(atol: Realtype) -> Self { AbsTolerance::Scalar(atol) } pub fn vector(atol: &[Realtype; SIZE]) -> Self { let atol = NVectorSerialHeapAllocated::new_from(atol); AbsTolerance::Vector(atol) } } /// An enum representing the choice between scalars or vectors absolute tolerances /// for sensitivities. pub enum SensiAbsTolerance<const SIZE: usize, const N_SENSI: usize> { Scalar([Realtype; N_SENSI]), Vector([NVectorSerialHeapAllocated<SIZE>; N_SENSI]), } impl<const SIZE: usize, const N_SENSI: usize> SensiAbsTolerance<SIZE, N_SENSI> { pub fn scalar(atol: [Realtype; N_SENSI]) -> Self { SensiAbsTolerance::Scalar(atol) } pub fn vector(atol: &[[Realtype; SIZE]; N_SENSI]) -> Self { SensiAbsTolerance::Vector( array_init::from_iter( atol.iter() .map(|arr| NVectorSerialHeapAllocated::new_from(arr)), ) .unwrap(), ) } } /// A short-hand for `std::result::Result<T, crate::Error>` pub type Result<T> = std::result::Result<T, Error>; fn check_non_null<T>(ptr: *mut T, func_id: &'static str) -> Result<NonNull<T>> { NonNull::new(ptr).ok_or(Error::NullPointerError { func_id }) } fn check_flag_is_succes(flag: c_int, func_id: &'static str) -> Result<()> { if flag == sundials_sys::CV_SUCCESS { Ok(()) } else { Err(Error::ErrorCode { flag, func_id }) } } #[repr(C)] struct CvodeMemoryBlock { _private: [u8; 0], } #[repr(transparent)] #[derive(Debug, Clone, Copy)] struct CvodeMemoryBlockNonNullPtr { ptr: NonNull<CvodeMemoryBlock>, } impl CvodeMemoryBlockNonNullPtr { fn new(ptr: NonNull<CvodeMemoryBlock>) -> Self { Self { ptr } } fn as_raw(self) -> *mut c_void { self.ptr.as_ptr() as *mut c_void } } impl From<NonNull<CvodeMemoryBlock>> for CvodeMemoryBlockNonNullPtr { fn from(x: NonNull<CvodeMemoryBlock>) -> Self { Self::new(x) } }