From 77506226b9a451f93a619dd2fd2b027d66e25a97 Mon Sep 17 00:00:00 2001 From: Dorival Pedroso Date: Fri, 8 Mar 2024 20:09:45 +1000 Subject: [PATCH] [ode] Add tests --- russell_ode/Cargo.toml | 5 +- russell_ode/README.md | 2 +- russell_ode/src/radau5.rs | 117 ++++++++++++++++++++++++++++++++++++- russell_ode/src/samples.rs | 18 ++++-- 4 files changed, 132 insertions(+), 10 deletions(-) diff --git a/russell_ode/Cargo.toml b/russell_ode/Cargo.toml index b0bd4350..d3227e92 100644 --- a/russell_ode/Cargo.toml +++ b/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" diff --git a/russell_ode/README.md b/russell_ode/README.md index 39b8ba8f..0e9c2e5b 100644 --- a/russell_ode/README.md +++ b/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)_ diff --git a/russell_ode/src/radau5.rs b/russell_ode/src/radau5.rs index 45a7df06..71ba11cc 100644 --- a/russell_ode/src/radau5.rs +++ b/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); + } + } } diff --git a/russell_ode/src/samples.rs b/russell_ode/src/samples.rs index cf3b0841..085cb23f 100644 --- a/russell_ode/src/samples.rs +++ b/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) }