Merge pull request #112 from cpmech/fix-ode-output

Fix ode output

russell_lab/build.rs

@@ -24,7 +24,6 @@ fn compile_blas() {
     println!("cargo:rustc-link-lib=m");
     println!("cargo:rustc-link-lib=dl");
     println!("cargo:rustc-link-lib=iomp5");
-    println!("cargo:rustc-cfg=use_intel_mkl");
 }
 
 // OpenBLAS

russell_lab/src/base/auxiliary_blas.rs

@@ -79,7 +79,7 @@ mod tests {
 
     #[test]
     fn using_intel_mkl_works() {
-        if cfg!(use_intel_mkl) {
+        if cfg!(feature = "intel_mkl") {
             assert!(using_intel_mkl());
         } else {
             assert!(!using_intel_mkl());

russell_ode/Cargo.toml

@@ -13,6 +13,8 @@ keywords = ["differential", "numerical", "solver"]
 
 [features]
 intel_mkl = ["russell_lab/intel_mkl", "russell_sparse/intel_mkl"]
+local_suitesparse = ["russell_sparse/local_suitesparse"]
+with_mumps = ["russell_sparse/with_mumps"]
 
 [dependencies]
 russell_lab = { path = "../russell_lab", version = "1.0.0" }

russell_ode/README.md

@@ -90,6 +90,8 @@ russell_ode = "*"
 The following (Rust) features are available:
 
 * `intel_mkl`: Use Intel MKL instead of OpenBLAS
+* `local_suitesparse`: Use a locally compiled version of SuiteSparse
+* `with_mumps`: Enable the MUMPS solver (locally compiled)
 
 Note that the [main README file](https://github.com/cpmech/russell) presents the steps to compile the required libraries according to each feature.
 
@@ -148,7 +150,7 @@ fn main() -> Result<(), StrError> {
     // solve from x = 0 to x = 1
     let x1 = 1.0;
     let mut args = 0;
-    solver.solve(&mut y, x, x1, None, None, &mut args)?;
+    solver.solve(&mut y, x, x1, None, &mut args)?;
     println!("y =\n{}", y);
 
     // check the results
@@ -280,7 +282,7 @@ fn main() -> Result<(), StrError> {
     // solve from x = 0 to x = 20
     let x1 = 20.0;
     let mut args = 0;
-    solver.solve(&mut y, x, x1, None, None, &mut args)?;
+    solver.solve(&mut y, x, x1, None, &mut args)?;
     println!("y =\n{}", y);
 
     // check the results

russell_ode/src/lib.rs

@@ -179,7 +179,7 @@
 //!     // solve from x = 0 to x = 2
 //!     let x1 = 2.0;
 //!     let mut args = 0;
-//!     solver.solve(&mut y, x, x1, None, None, &mut args)?;
+//!     solver.solve(&mut y, x, x1, None, &mut args)?;
 //!     println!("y =\n{}", y);
 //!     Ok(())
 //! }

russell_ode/src/ode_solver.rs

@@ -99,7 +99,7 @@ use russell_sparse::CooMatrix;
 ///     // solve from x = 0 to x = 1
 ///     let x1 = 1.0;
 ///     let mut args = 0;
-///     solver.solve(&mut y, x, x1, None, None, &mut args)?;
+///     solver.solve(&mut y, x, x1, None, &mut args)?;
 ///
 ///     // check the results
 ///     let y_ana = Vector::from(&[f64::exp(x1) - x1 - 1.0]);
@@ -432,7 +432,7 @@ impl<'a, A> OdeSolver<'a, A> {
 #[cfg(test)]
 mod tests {
     use super::OdeSolver;
-    use crate::{no_jacobian, HasJacobian, NoArgs, OutCallback, OutCount, OutData};
+    use crate::{no_jacobian, HasJacobian, NoArgs, OutCount, OutData, Stats};
     use crate::{Method, Params, Samples, System};
     use russell_lab::{approx_eq, array_approx_eq, vec_approx_eq, Vector};
     use russell_sparse::Genie;
@@ -628,7 +628,7 @@ mod tests {
         }
 
         // define the callback function
-        let cb: OutCallback<NoArgs> = |_stats, _h, _x, _y, _args| Err("unreachable");
+        let cb = |_stats: &Stats, _h: f64, _x: f64, _y: &Vector, _args: &mut NoArgs| Err("unreachable");
         assert_eq!(cb(&solver.stats(), 0.0, 0.0, &y0, &mut args).err(), Some("unreachable"));
 
         // run again without step output
@@ -736,7 +736,7 @@ mod tests {
         }
 
         // define the callback function
-        let cb: OutCallback<NoArgs> = |_stats, _h, _x, _y, _args| Err("unreachable");
+        let cb = |_stats: &Stats, _h: f64, _x: f64, _y: &Vector, _args: &mut NoArgs| Err("unreachable");
         assert_eq!(cb(&solver.stats(), 0.0, 0.0, &y0, &mut args).err(), Some("unreachable"));
 
         // run again without dense output

russell_ode/src/output.rs

@@ -8,18 +8,6 @@ use std::io::BufReader;
 use std::marker::PhantomData;
 use std::path::Path;
 
-/// Defines a function to compute y(x) (e.g, the analytical solution)
-///
-/// Use `|y, x, args|` or `|y: &mut Vector, x: f64, args, &mut A|`
-pub type YxFunction<A> = fn(&mut Vector, f64, &mut A);
-
-/// Defines a callback function to be executed during the output of results (accepted step or dense output)
-///
-/// Use `|stats, h, x, y, args|` or `|stats: &Stats, h: f64, x: f64, y: &Vector, args: &mut A|`
-///
-/// The function may return `true` to stop the computations
-pub type OutCallback<A> = fn(&Stats, f64, f64, &Vector, &mut A) -> Result<bool, StrError>;
-
 /// Holds the data generated at an accepted step or during the dense output
 #[derive(Clone, Debug, Deserialize)]
 pub struct OutData {
@@ -95,7 +83,7 @@ pub struct Output<'a, A> {
     dense_last_x: f64,
 
     /// Holds a callback function for the dense output
-    dense_callback: Option<OutCallback<A>>,
+    dense_callback: Option<Box<dyn Fn(&Stats, f64, f64, &Vector, &mut A) -> Result<bool, StrError> + 'a>>,
 
     /// Save the results to a file (dense)
     dense_file_key: Option<String>,
@@ -136,7 +124,7 @@ pub struct Output<'a, A> {
     y_aux: Vector,
 
     /// Holds the y(x) function (e.g., to compute the correct/analytical solution)
-    yx_function: Option<YxFunction<A>>,
+    yx_function: Option<Box<dyn Fn(&mut Vector, f64, &mut A) + 'a>>,
 
     /// Handles the generic argument
     phantom: PhantomData<fn() -> A>,
@@ -302,14 +290,14 @@ impl<'a, A> Output<'a, A> {
         &mut self,
         enable: bool,
         h_out: f64,
-        callback: OutCallback<A>,
+        callback: impl Fn(&Stats, f64, f64, &Vector, &mut A) -> Result<bool, StrError> + 'a,
     ) -> Result<&mut Self, StrError> {
         if h_out <= f64::EPSILON {
             return Err("h_out must be > EPSILON");
         }
         self.dense_h_out = h_out;
         if enable {
-            self.dense_callback = Some(callback);
+            self.dense_callback = Some(Box::new(callback));
         } else {
             self.dense_callback = None;
         }
@@ -381,8 +369,10 @@ impl<'a, A> Output<'a, A> {
     }
 
     /// Sets the function to compute the correct/reference results y(x)
-    pub fn set_yx_correct(&mut self, y_fn_x: fn(&mut Vector, f64, &mut A)) -> &mut Self {
-        self.yx_function = Some(y_fn_x);
+    ///
+    /// Use `|y, x, args|` or `|y: &mut Vector, x: f64, args, &mut A|`
+    pub fn set_yx_correct(&mut self, y_fn_x: impl Fn(&mut Vector, f64, &mut A) + 'a) -> &mut Self {
+        self.yx_function = Some(Box::new(y_fn_x));
         self
     }
 
@@ -423,7 +413,7 @@ impl<'a, A> Output<'a, A> {
         //
         // step output: callback
         if let Some(cb) = self.step_callback.as_ref() {
-            let stop = (cb)(&work.stats, h, x, y, args)?;
+            let stop = cb(&work.stats, h, x, y, args)?;
             if stop {
                 return Ok(stop);
             }
@@ -462,7 +452,7 @@ impl<'a, A> Output<'a, A> {
                 self.dense_last_x = x;
 
                 // first dense output: callback
-                if let Some(cb) = self.dense_callback {
+                if let Some(cb) = self.dense_callback.as_ref() {
                     let stop = cb(&work.stats, h, x, y, args)?;
                     if stop {
                         return Ok(stop);
@@ -502,7 +492,7 @@ impl<'a, A> Output<'a, A> {
                     solver.dense_output(y_out, x_out, x, y, h);
 
                     // subsequent dense output: callback
-                    if let Some(cb) = self.dense_callback {
+                    if let Some(cb) = self.dense_callback.as_ref() {
                         let stop = cb(&work.stats, h, x_out, y_out, args)?;
                         if stop {
                             return Ok(stop);
@@ -561,7 +551,7 @@ impl<'a, A> Output<'a, A> {
         }
 
         // dense output: callback
-        if let Some(cb) = self.dense_callback {
+        if let Some(cb) = self.dense_callback.as_ref() {
             cb(&work.stats, h, x, y, args)?;
         }
 

russell_ode/src/samples.rs

@@ -1,5 +1,5 @@
 use crate::StrError;
-use crate::{HasJacobian, NoArgs, PdeDiscreteLaplacian2d, Side, System, YxFunction};
+use crate::{HasJacobian, NoArgs, PdeDiscreteLaplacian2d, Side, System};
 use russell_lab::math::PI;
 use russell_lab::Vector;
 use russell_sparse::{CooMatrix, Genie, Sym};
@@ -53,7 +53,7 @@ impl Samples {
         f64,
         Vector,
         NoArgs,
-        YxFunction<NoArgs>,
+        impl Fn(&mut Vector, f64, &mut NoArgs),
     ) {
         // system
         let ndim = 1;
@@ -169,7 +169,7 @@ impl Samples {
         f64,
         Vector,
         NoArgs,
-        YxFunction<NoArgs>,
+        impl Fn(&mut Vector, f64, &mut NoArgs),
     ) {
         // selected symmetric option (for both Mass and Jacobian matrices)
         let sym = genie.symmetry(symmetric);
@@ -820,7 +820,7 @@ impl Samples {
         f64,
         Vector,
         NoArgs,
-        YxFunction<NoArgs>,
+        impl Fn(&mut Vector, f64, &mut NoArgs),
     ) {
         // constants
         const L: f64 = -50.0; // lambda
@@ -1225,7 +1225,7 @@ impl Samples {
         f64,
         Vector,
         NoArgs,
-        YxFunction<NoArgs>,
+        impl Fn(&mut Vector, f64, &mut NoArgs),
     ) {
         // system
         let ndim = 1;
@@ -1306,7 +1306,7 @@ impl Samples {
         f64,
         Vector,
         NoArgs,
-        YxFunction<NoArgs>,
+        impl Fn(&mut Vector, f64, &mut NoArgs),
     ) {
         // system
         let ndim = 2;

russell_ode/tests/test_radau5_robertson_small_h.rs

@@ -14,12 +14,12 @@ fn test_radau5_robertson_small_h() {
     params.step.h_ini = 1e-6;
     params.debug = true;
 
-    // allocate the solver
-    let mut solver = OdeSolver::new(params, &system).unwrap();
-
     // this will cause h to become too small
     params.set_tolerances(1e-2, 1e-2, None).unwrap();
 
+    // allocate the solver
+    let mut solver = OdeSolver::new(params, &system).unwrap();
+
     // solve the ODE system
     let res = solver.solve(&mut y0, x0, x1, None, &mut args);
     assert_eq!(res.err(), Some("the stepsize becomes too small"));

russell_sparse/Cargo.toml

@@ -12,9 +12,9 @@ categories = ["mathematics", "science"]
 keywords = ["matrix", "sparse", "solver"]
 
 [features]
+intel_mkl = ["russell_lab/intel_mkl"]
 local_suitesparse = []
 with_mumps = []
-intel_mkl = ["russell_lab/intel_mkl"]
 
 [dependencies]
 num-traits = "0.2"

russell_sparse/README.md

@@ -63,9 +63,9 @@ russell_sparse = "*"
 
 The following (Rust) features are available:
 
+* `intel_mkl`: Use Intel MKL instead of OpenBLAS
 * `local_suitesparse`: Use a locally compiled version of SuiteSparse
 * `with_mumps`: Enable the MUMPS solver (locally compiled)
-* `intel_mkl`: Use Intel MKL instead of OpenBLAS
 
 Note that the [main README file](https://github.com/cpmech/russell) presents the steps to compile the required libraries according to each feature.