[ode] Add tests

russell_ode/Cargo.toml

@@ -3,7 +3,7 @@ name = "russell_ode"
 version = "0.1.0"
 edition = "2021"
 license = "MIT"
-description = "Solvers for Ordinary Differential Equations"
+description = "Solvers for Ordinary Differential Equations and Differential Algebraic Systems"
 homepage = "https://github.com/cpmech/russell"
 repository = "https://github.com/cpmech/russell"
 documentation = "https://docs.rs/russell_ode"
@@ -15,3 +15,6 @@ keywords = ["differential equations", "numerical methods", "solver"]
 russell_lab = { path = "../russell_lab", version = "0.8" }
 russell_sparse = { path = "../russell_sparse", version = "0.8" }
 num-complex = { version = "0.4", features = ["serde"] }
+
+[dev-dependencies]
+serial_test = "3.0"

russell_ode/README.md

@@ -1,4 +1,4 @@
-# Russell ODE - Solvers for Ordinary Differential Equations
+# Russell ODE - Solvers for Ordinary Differential Equations and Differential Algebraic Systems
 
 _This crate is part of [Russell - Rust Scientific Library](https://github.com/cpmech/russell)_
 

russell_ode/src/radau5.rs

@@ -701,7 +701,12 @@ const TI: [[f64; 3]; 3] = [
 mod tests {
     use super::Radau5;
     use crate::{Method, OdeSolverTrait, Params, Samples, Workspace};
-    use russell_lab::{format_fortran, Vector};
+    use russell_lab::{format_fortran, format_scientific, Vector};
+    use russell_sparse::Genie;
+    use serial_test::serial;
+
+    // IMPORTANT:
+    // Since MUMPS is not thread-safe, we need to use serial_test::serial
 
     #[test]
     fn radau5_works() {
@@ -715,7 +720,7 @@ mod tests {
         // allocate structs
         let params = Params::new(Method::Radau5);
         let mut solver = Radau5::new(params, &system);
-        let mut work = Workspace::new(Method::FwEuler);
+        let mut work = Workspace::new(Method::Radau5);
 
         // message
         println!("{:>4}{:>23}{:>23}", "step", "err_y0", "err_y1");
@@ -785,7 +790,7 @@ mod tests {
         let mut params = Params::new(Method::Radau5);
         params.newton.use_numerical_jacobian = true;
         let mut solver = Radau5::new(params, &system);
-        let mut work = Workspace::new(Method::FwEuler);
+        let mut work = Workspace::new(Method::Radau5);
 
         // message
         println!("{:>4}{:>23}{:>23}", "step", "err_y0", "err_y1");
@@ -841,4 +846,110 @@ mod tests {
         // check number of function evaluations
         assert_eq!(work.bench.n_function, n_fcn_correct);
     }
+
+    #[test]
+    fn radau5_works_mass_matrix() {
+        // problem
+        let (system, data, mut args) = Samples::simple_system_with_mass_matrix(false);
+        let yfx = data.y_analytical.unwrap();
+        let ndim = system.ndim;
+
+        // allocate structs
+        let params = Params::new(Method::Radau5);
+        let mut solver = Radau5::new(params, &system);
+        let mut work = Workspace::new(Method::Radau5);
+
+        // message
+        println!("{:>4}{:>10}{:>10}{:>10}", "step", "err_y0", "err_y1", "err_y2");
+
+        // numerical approximation
+        let h = 0.1;
+        let mut x = data.x0;
+        let mut y = data.y0.clone();
+        let mut y_ana = Vector::new(ndim);
+        for n in 0..5 {
+            // call step
+            solver.step(&mut work, x, &y, h, &mut args).unwrap();
+
+            // important: update n_accepted (must precede `accept`)
+            work.bench.n_accepted += 1;
+
+            // call accept
+            solver.accept(&mut work, &mut x, &mut y, h, &mut args).unwrap();
+
+            // important: save previous stepsize and relative error (must succeed `accept`)
+            work.h_prev = h;
+            work.rel_error_prev = f64::max(params.step.rel_error_prev_min, work.rel_error);
+
+            // check the results
+            yfx(&mut y_ana, x);
+            let err_y0 = f64::abs(y[0] - y_ana[0]);
+            let err_y1 = f64::abs(y[1] - y_ana[1]);
+            let err_y2 = f64::abs(y[2] - y_ana[2]);
+            println!(
+                "{:>4}{}{}{}",
+                n,
+                format_scientific(err_y0, 10, 2),
+                format_scientific(err_y1, 10, 2),
+                format_scientific(err_y2, 10, 2)
+            );
+            assert!(err_y0 < 1e-9);
+            assert!(err_y1 < 1e-9);
+            assert!(err_y2 < 1e-8);
+        }
+    }
+
+    #[test]
+    #[serial]
+    fn radau5_works_mass_matrix_symmetric_mumps() {
+        // problem
+        let (system, data, mut args) = Samples::simple_system_with_mass_matrix(true);
+        let yfx = data.y_analytical.unwrap();
+        let ndim = system.ndim;
+
+        // allocate structs
+        let mut params = Params::new(Method::Radau5);
+        params.newton.genie = Genie::Mumps;
+        let mut solver = Radau5::new(params, &system);
+        let mut work = Workspace::new(Method::Radau5);
+
+        // message
+        println!("{:>4}{:>10}{:>10}{:>10}", "step", "err_y0", "err_y1", "err_y2");
+
+        // numerical approximation
+        let h = 0.1;
+        let mut x = data.x0;
+        let mut y = data.y0.clone();
+        let mut y_ana = Vector::new(ndim);
+        for n in 0..5 {
+            // call step
+            solver.step(&mut work, x, &y, h, &mut args).unwrap();
+
+            // important: update n_accepted (must precede `accept`)
+            work.bench.n_accepted += 1;
+
+            // call accept
+            solver.accept(&mut work, &mut x, &mut y, h, &mut args).unwrap();
+
+            // important: save previous stepsize and relative error (must succeed `accept`)
+            work.h_prev = h;
+            work.rel_error_prev = f64::max(params.step.rel_error_prev_min, work.rel_error);
+
+            // check the results
+            yfx(&mut y_ana, x);
+            let err_y0 = f64::abs(y[0] - y_ana[0]);
+            let err_y1 = f64::abs(y[1] - y_ana[1]);
+            let err_y2 = f64::abs(y[2] - y_ana[2]);
+            println!(
+                "{:>4}{}{}{}",
+                n,
+                format_scientific(err_y0, 10, 2),
+                format_scientific(err_y1, 10, 2),
+                format_scientific(err_y2, 10, 2)
+            );
+            assert!(err_y0 < 1e-9);
+            assert!(err_y1 < 1e-9);
+            assert!(err_y2 < 1e-8);
+        }
+    }
 }

russell_ode/src/samples.rs

@@ -176,7 +176,7 @@ impl Samples {
 
         // ODE system
         let ndim = 3;
-        let jac_nnz = 4;
+        let jac_nnz = if lower_triangle { 3 } else { 4 };
         let mut system = System::new(
             ndim,
             move |f: &mut Vector, x: f64, y: &Vector, _args: &mut NoArgs| {
@@ -188,7 +188,9 @@ impl Samples {
             move |jj: &mut CooMatrix, _x: f64, _y: &Vector, m: f64, _args: &mut NoArgs| {
                 jj.reset();
                 jj.put(0, 0, m * (-1.0)).unwrap();
-                jj.put(0, 1, m * (1.0)).unwrap();
+                if !lower_triangle {
+                    jj.put(0, 1, m * (1.0)).unwrap();
+                }
                 jj.put(1, 0, m * (1.0)).unwrap();
                 jj.put(1, 1, m * (1.0)).unwrap();
                 Ok(())
@@ -199,10 +201,12 @@ impl Samples {
         );
 
         // mass matrix
-        let mass_nnz = 5;
+        let mass_nnz = if lower_triangle { 4 } else { 5 };
         system.init_mass_matrix(mass_nnz).unwrap();
         system.mass_put(0, 0, 1.0).unwrap();
-        system.mass_put(0, 1, 1.0).unwrap();
+        if !lower_triangle {
+            system.mass_put(0, 1, 1.0).unwrap();
+        }
         system.mass_put(1, 0, 1.0).unwrap();
         system.mass_put(1, 1, -1.0).unwrap();
         system.mass_put(2, 2, 1.0).unwrap();
@@ -213,7 +217,11 @@ impl Samples {
             y0,
             x1,
             h_equal: None,
-            y_analytical: None,
+            y_analytical: Some(|y, x| {
+                y[0] = f64::cos(x);
+                y[1] = -f64::sin(x);
+                y[2] = f64::ln(1.0 + x);
+            }),
         };
         (system, data, 0)
     }