From 70bf0df63ca1998f00230e08cc422b7a55acb13a Mon Sep 17 00:00:00 2001 From: Dorival Pedroso Date: Wed, 22 May 2024 14:23:55 +1000 Subject: [PATCH 1/2] Replace the mass matrix in System by a callback function. Store mass matrix in Radau5. --- russell_ode/README.md | 13 +- .../examples/simple_system_with_mass.rs | 13 +- russell_ode/src/lib.rs | 13 +- russell_ode/src/ode_solver.rs | 2 +- russell_ode/src/radau5.rs | 22 ++-- russell_ode/src/samples.rs | 46 ++++--- russell_ode/src/system.rs | 116 +++++++++++------- russell_ode/zscripts/memcheck.bash | 2 +- 8 files changed, 134 insertions(+), 93 deletions(-) diff --git a/russell_ode/README.md b/russell_ode/README.md index 26a3ca2f..4642dd19 100644 --- a/russell_ode/README.md +++ b/russell_ode/README.md @@ -259,12 +259,13 @@ fn main() -> Result<(), StrError> { // mass matrix let mass_nnz = 5; - system.init_mass_matrix(mass_nnz, symmetric)?; - system.mass_put(0, 0, 1.0)?; - system.mass_put(0, 1, 1.0)?; - system.mass_put(1, 0, 1.0)?; - system.mass_put(1, 1, -1.0)?; - system.mass_put(2, 2, 1.0)?; + system.set_mass(Some(mass_nnz), symmetric, |mm: &mut CooMatrix| { + mm.put(0, 0, 1.0).unwrap(); + mm.put(0, 1, 1.0).unwrap(); + mm.put(1, 0, 1.0).unwrap(); + mm.put(1, 1, -1.0).unwrap(); + mm.put(2, 2, 1.0).unwrap(); + })?; // solver let params = Params::new(Method::Radau5); diff --git a/russell_ode/examples/simple_system_with_mass.rs b/russell_ode/examples/simple_system_with_mass.rs index fdc4dd18..ee2cca3d 100644 --- a/russell_ode/examples/simple_system_with_mass.rs +++ b/russell_ode/examples/simple_system_with_mass.rs @@ -30,12 +30,13 @@ fn main() -> Result<(), StrError> { // mass matrix let mass_nnz = 5; - system.init_mass_matrix(mass_nnz, symmetric)?; - system.mass_put(0, 0, 1.0)?; - system.mass_put(0, 1, 1.0)?; - system.mass_put(1, 0, 1.0)?; - system.mass_put(1, 1, -1.0)?; - system.mass_put(2, 2, 1.0)?; + system.set_mass(Some(mass_nnz), symmetric, |mm: &mut CooMatrix| { + mm.put(0, 0, 1.0).unwrap(); + mm.put(0, 1, 1.0).unwrap(); + mm.put(1, 0, 1.0).unwrap(); + mm.put(1, 1, -1.0).unwrap(); + mm.put(2, 2, 1.0).unwrap(); + })?; // solver let params = Params::new(Method::Radau5); diff --git a/russell_ode/src/lib.rs b/russell_ode/src/lib.rs index 19bd0693..f3e1047b 100644 --- a/russell_ode/src/lib.rs +++ b/russell_ode/src/lib.rs @@ -162,12 +162,13 @@ //! //! // mass matrix //! let mass_nnz = 5; // number of non-zero values in the mass matrix -//! system.init_mass_matrix(mass_nnz, symmetric)?; -//! system.mass_put(0, 0, 1.0)?; -//! system.mass_put(0, 1, 1.0)?; -//! system.mass_put(1, 0, 1.0)?; -//! system.mass_put(1, 1, -1.0)?; -//! system.mass_put(2, 2, 1.0)?; +//! system.set_mass(Some(mass_nnz), symmetric, |mm: &mut CooMatrix| { +//! mm.put(0, 0, 1.0).unwrap(); +//! mm.put(0, 1, 1.0).unwrap(); +//! mm.put(1, 0, 1.0).unwrap(); +//! mm.put(1, 1, -1.0).unwrap(); +//! mm.put(2, 2, 1.0).unwrap(); +//! })?; //! //! // solver //! let params = Params::new(Method::Radau5); diff --git a/russell_ode/src/ode_solver.rs b/russell_ode/src/ode_solver.rs index 0275a0f7..963a93ad 100644 --- a/russell_ode/src/ode_solver.rs +++ b/russell_ode/src/ode_solver.rs @@ -139,7 +139,7 @@ impl<'a, A> OdeSolver<'a, A> { where A: 'a, { - if system.mass_matrix.is_some() && params.method != Method::Radau5 { + if system.calc_mass.is_some() && params.method != Method::Radau5 { return Err("the method must be Radau5 for systems with a mass matrix"); } params.validate()?; diff --git a/russell_ode/src/radau5.rs b/russell_ode/src/radau5.rs index 54dbda1c..8663ac97 100644 --- a/russell_ode/src/radau5.rs +++ b/russell_ode/src/radau5.rs @@ -2,7 +2,7 @@ use crate::StrError; use crate::{OdeSolverTrait, Params, System, Workspace}; use russell_lab::math::SQRT_6; use russell_lab::{complex_vec_zip, cpx, format_fortran, vec_copy, Complex64, ComplexVector, Vector}; -use russell_sparse::{numerical_jacobian, ComplexCscMatrix, CscMatrix}; +use russell_sparse::{numerical_jacobian, ComplexCscMatrix, CooMatrix, CscMatrix}; use russell_sparse::{ComplexLinSolver, ComplexSparseMatrix, Genie, LinSolver, SparseMatrix}; use std::thread; @@ -31,6 +31,9 @@ pub(crate) struct Radau5<'a, A> { /// ODE system system: &'a System<'a, A>, + /// Holds the mass matrix + mass: Option, + /// Holds the Jacobian matrix. J = df/dy jj: SparseMatrix, @@ -119,9 +122,13 @@ impl<'a, A> Radau5<'a, A> { /// Allocates a new instance pub fn new(params: Params, system: &'a System<'a, A>) -> Self { let ndim = system.ndim; - let mass_nnz = match system.mass_matrix.as_ref() { - Some(mass) => mass.get_info().2, - None => ndim, + let (mass, mass_nnz) = match system.calc_mass.as_ref() { + Some(calc) => { + let mut mm = CooMatrix::new(ndim, ndim, system.mass_nnz, system.symmetric).unwrap(); + (calc)(&mut mm); + (Some(mm), system.mass_nnz) + } + None => (None, ndim), // ndim => diagonal }; let jac_nnz = if params.newton.use_numerical_jacobian { if system.symmetric.triangular() { @@ -137,6 +144,7 @@ impl<'a, A> Radau5<'a, A> { Radau5 { params, system, + mass, jj: SparseMatrix::new_coo(ndim, ndim, jac_nnz, system.symmetric).unwrap(), kk_real: SparseMatrix::new_coo(ndim, ndim, nnz, system.symmetric).unwrap(), kk_comp: ComplexSparseMatrix::new_coo(ndim, ndim, nnz, system.symmetric).unwrap(), @@ -214,7 +222,7 @@ impl<'a, A> Radau5<'a, A> { let gamma = GAMMA / h; kk_real.assign(-1.0, jj).unwrap(); // K_real = -J kk_comp.assign_real(-1.0, 0.0, jj).unwrap(); // K_comp = -J - match self.system.mass_matrix.as_ref() { + match self.mass.as_ref() { Some(mass) => { kk_real.augment(gamma, mass).unwrap(); // K_real += γ M kk_comp.augment_real(alpha, beta, mass).unwrap(); // K_comp += (α + βi) M @@ -418,7 +426,7 @@ impl<'a, A> OdeSolverTrait for Radau5<'a, A> { (self.system.function)(&mut self.k2, u2, &self.v2, args)?; // compute the right-hand side vectors - let (l0, l1, l2) = match self.system.mass_matrix.as_ref() { + let (l0, l1, l2) = match self.mass.as_ref() { Some(mass) => { mass.mat_vec_mul(&mut self.dw0, 1.0, &self.w0).unwrap(); // dw0 := M ⋅ w0 mass.mat_vec_mul(&mut self.dw1, 1.0, &self.w1).unwrap(); // dw1 := M ⋅ w1 @@ -538,7 +546,7 @@ impl<'a, A> OdeSolverTrait for Radau5<'a, A> { let err = &mut self.dw0; // error variable // compute ez, mez and rhs - match self.system.mass_matrix.as_ref() { + match self.mass.as_ref() { Some(mass) => { for m in 0..ndim { ez[m] = E0 * self.z0[m] + E1 * self.z1[m] + E2 * self.z2[m]; diff --git a/russell_ode/src/samples.rs b/russell_ode/src/samples.rs index ca8b7da5..9a8b1fd0 100644 --- a/russell_ode/src/samples.rs +++ b/russell_ode/src/samples.rs @@ -192,14 +192,17 @@ impl Samples { // mass matrix let mass_nnz = if triangular { 4 } else { 5 }; - system.init_mass_matrix(mass_nnz, sym).unwrap(); - system.mass_put(0, 0, 1.0).unwrap(); - if !triangular { - 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(); + system + .set_mass(Some(mass_nnz), sym, move |mm: &mut CooMatrix| { + mm.put(0, 0, 1.0).unwrap(); + if !triangular { + mm.put(0, 1, 1.0).unwrap(); + } + mm.put(1, 0, 1.0).unwrap(); + mm.put(1, 1, -1.0).unwrap(); + mm.put(2, 2, 1.0).unwrap(); + }) + .unwrap(); // initial values let x0 = 0.0; @@ -1088,16 +1091,19 @@ impl Samples { const C2: f64 = 2e-6; const C3: f64 = 3e-6; let mass_nnz = 9; - system.init_mass_matrix(mass_nnz, symmetric).unwrap(); - system.mass_put(0, 0, -C1).unwrap(); - system.mass_put(0, 1, C1).unwrap(); - system.mass_put(1, 0, C1).unwrap(); - system.mass_put(1, 1, -C1).unwrap(); - system.mass_put(2, 2, -C2).unwrap(); - system.mass_put(3, 3, -C3).unwrap(); - system.mass_put(3, 4, C3).unwrap(); - system.mass_put(4, 3, C3).unwrap(); - system.mass_put(4, 4, -C3).unwrap(); + system + .set_mass(Some(mass_nnz), symmetric, |mm: &mut CooMatrix| { + mm.put(0, 0, -C1).unwrap(); + mm.put(0, 1, C1).unwrap(); + mm.put(1, 0, C1).unwrap(); + mm.put(1, 1, -C1).unwrap(); + mm.put(2, 2, -C2).unwrap(); + mm.put(3, 3, -C3).unwrap(); + mm.put(3, 4, C3).unwrap(); + mm.put(4, 3, C3).unwrap(); + mm.put(4, 4, -C3).unwrap(); + }) + .unwrap(); // initial values let x0 = 0.0; @@ -1531,7 +1537,9 @@ mod tests { mat_approx_eq(&ana, &num, 1e-13); // check the mass matrix - let mass = system.mass_matrix.unwrap(); + let mut mass = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.symmetric).unwrap(); + let calc_mass = system.calc_mass.as_ref().unwrap(); + (calc_mass)(&mut mass); println!("{}", mass.as_dense()); let ndim = system.ndim; let nnz_mass = 5 + 4; diff --git a/russell_ode/src/system.rs b/russell_ode/src/system.rs index 14354ed0..6fa93e59 100644 --- a/russell_ode/src/system.rs +++ b/russell_ode/src/system.rs @@ -71,9 +71,15 @@ pub struct System<'a, A> { /// Jacobian function pub(crate) jacobian: Option Result<(), StrError> + 'a>>, + /// Calc mass matrix function + pub(crate) calc_mass: Option>, + /// Number of non-zeros in the Jacobian matrix pub(crate) jac_nnz: usize, + /// Number of non-zeros in the mass matrix + pub(crate) mass_nnz: usize, + /// Symmetric type of the Jacobian matrix (for error checking; to make sure it is equal to sym_mass) sym_jac: Option, @@ -83,9 +89,6 @@ pub struct System<'a, A> { /// Symmetric type of the Jacobian and mass matrices pub(crate) symmetric: Sym, - /// Holds the mass matrix - pub(crate) mass_matrix: Option, - /// Handle generic argument phantom: PhantomData A>, } @@ -156,11 +159,12 @@ impl<'a, A> System<'a, A> { ndim, function: Box::new(function), jacobian: None, + calc_mass: None, jac_nnz: ndim * ndim, + mass_nnz: 0, sym_jac: None, sym_mass: None, symmetric: Sym::No, - mass_matrix: None, phantom: PhantomData, } } @@ -171,14 +175,14 @@ impl<'a, A> System<'a, A> { /// /// # Input /// - /// * `jac_nnz` -- the number of non-zeros in the Jacobian; use None to indicate a dense matrix with: + /// * `nnz` -- the number of non-zeros in the Jacobian; use None to indicate a dense matrix with: /// * `nnz = (ndim + ndim²) / 2` if triangular /// * `nnz = ndim²` otherwise /// * `symmetric` -- specifies the symmetric type of the Jacobian and **mass** matrices /// * `callback` -- the function to calculate the Jacobian matrix pub fn set_jacobian( &mut self, - jac_nnz: Option, + nnz: Option, symmetric: Sym, callback: impl Fn(&mut CooMatrix, f64, f64, &Vector, &mut A) -> Result<(), StrError> + 'a, ) -> Result<(), StrError> { @@ -187,8 +191,8 @@ impl<'a, A> System<'a, A> { return Err("the Jacobian matrix must have the same symmetric type as the mass matrix"); } } - self.jac_nnz = if let Some(nnz) = jac_nnz { - nnz + self.jac_nnz = if let Some(value) = nnz { + value } else { if symmetric.triangular() { (self.ndim + self.ndim * self.ndim) / 2 @@ -202,45 +206,41 @@ impl<'a, A> System<'a, A> { Ok(()) } - /// Initializes and enables the mass matrix - /// - /// **Note:** Even if the (analytical) Jacobian function is not configured, - /// a numerical Jacobian matrix may be computed (see [crate::Params] and [crate::ParamsNewton]). + /// Sets a function to calculate the constant mass matrix /// /// # Input /// - /// * `max_nnz` -- max number of non-zero values + /// * `nnz` -- the number of non-zeros in the mass matrix; use None to indicate a dense matrix with: + /// * `nnz = (ndim + ndim²) / 2` if triangular + /// * `nnz = ndim²` otherwise /// * `symmetric` -- specifies the symmetric type for the mass and **Jacobian** matrices - /// - /// Use [System::mass_put] to "put" elements into the mass matrix. - pub fn init_mass_matrix(&mut self, max_nnz: usize, symmetric: Sym) -> Result<(), StrError> { + /// * `callback` -- the function to calculate the mass matrix (will be called just once) + pub fn set_mass( + &mut self, + nnz: Option, + symmetric: Sym, + callback: impl Fn(&mut CooMatrix) + 'a, + ) -> Result<(), StrError> { if let Some(sym) = self.sym_jac { if symmetric != sym { return Err("the mass matrix must have the same symmetric type as the Jacobian matrix"); } } + self.mass_nnz = if let Some(value) = nnz { + value + } else { + if symmetric.triangular() { + (self.ndim + self.ndim * self.ndim) / 2 + } else { + self.ndim * self.ndim + } + }; self.sym_mass = Some(symmetric); self.symmetric = symmetric; - self.mass_matrix = Some(CooMatrix::new(self.ndim, self.ndim, max_nnz, self.symmetric).unwrap()); + self.calc_mass = Some(Box::new(callback)); Ok(()) } - /// Puts a new element in the mass matrix (duplicates allowed) - /// - /// See also [russell_sparse::CooMatrix::put]. - /// - /// # Input - /// - /// * `i` -- row index (indices start at zero; zero-based) - /// * `j` -- column index (indices start at zero; zero-based) - /// * `value` -- the value M(i,j) - pub fn mass_put(&mut self, i: usize, j: usize, value: f64) -> Result<(), StrError> { - match self.mass_matrix.as_mut() { - Some(mass) => mass.put(i, j, value), - None => Err("mass matrix has not been initialized/enabled"), - } - } - /// Returns the dimension of the ODE system pub fn get_ndim(&self) -> usize { self.ndim @@ -250,6 +250,11 @@ impl<'a, A> System<'a, A> { pub fn get_jac_nnz(&self) -> usize { self.jac_nnz } + + /// Returns the number of non-zero values in the mass matrix + pub fn get_mass_nnz(&self) -> usize { + self.mass_nnz + } } //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -272,26 +277,25 @@ mod tests { let y = Vector::new(1); let mut args = 0; (system.function)(&mut f, x, &y, &mut args).unwrap(); - assert_eq!( - system.mass_put(0, 0, 1.0).err(), - Some("mass matrix has not been initialized/enabled") - ); - let cb = |_: &mut CooMatrix, _: f64, _: f64, _: &Vector, _: &mut NoArgs| Ok(()); + let jac_cb = |_: &mut CooMatrix, _: f64, _: f64, _: &Vector, _: &mut NoArgs| Ok(()); + let mas_cb = |_: &mut CooMatrix| (); let mut jj = CooMatrix::new(1, 1, 1, Sym::YesLower).unwrap(); + let mut mm = CooMatrix::new(1, 1, 1, Sym::YesLower).unwrap(); let y = Vector::new(1); - (cb)(&mut jj, 0.0, 0.0, &y, &mut 0).unwrap(); - system.set_jacobian(None, Sym::YesLower, cb).unwrap(); + (jac_cb)(&mut jj, 0.0, 0.0, &y, &mut 0).unwrap(); + (mas_cb)(&mut mm); + system.set_jacobian(None, Sym::YesLower, jac_cb).unwrap(); assert_eq!( - system.init_mass_matrix(1, Sym::YesUpper).err(), + system.set_mass(None, Sym::YesUpper, mas_cb).err(), Some("the mass matrix must have the same symmetric type as the Jacobian matrix") ); system.sym_jac = None; - system.init_mass_matrix(1, Sym::YesLower).unwrap(); + system.set_mass(None, Sym::YesLower, mas_cb).unwrap(); assert_eq!( - system.set_jacobian(None, Sym::YesUpper, cb).err(), + system.set_jacobian(None, Sym::YesUpper, jac_cb).err(), Some("the Jacobian matrix must have the same symmetric type as the mass matrix") ); - system.set_jacobian(None, Sym::YesLower, cb).unwrap(); // ok + system.set_jacobian(None, Sym::YesLower, jac_cb).unwrap(); // ok } #[test] @@ -361,9 +365,6 @@ mod tests { // analytical_solution: // y[0] = f64::cos(x * x / 2.0) - 2.0 * f64::sin(x * x / 2.0); // y[1] = 2.0 * f64::cos(x * x / 2.0) + f64::sin(x * x / 2.0); - system.init_mass_matrix(2, symmetric).unwrap(); // diagonal mass matrix => OK, but not needed - system.mass_put(0, 0, 1.0).unwrap(); - system.mass_put(1, 1, 1.0).unwrap(); // call system function let x = 0.0; let y = Vector::new(2); @@ -381,4 +382,25 @@ mod tests { assert_eq!(args.more_data_goes_here_fn, true); assert_eq!(args.more_data_goes_here_jj, true); } + + #[test] + fn ode_system_set_mass_works() { + let mut system = System::new(2, |f, _, _, _: &mut NoArgs| { + f[0] = 1.0; + f[1] = 1.0; + Ok(()) + }); + let mut f = Vector::new(2); + let x = 0.0; + let y = Vector::new(2); + let mut args = 0; + (system.function)(&mut f, x, &y, &mut args).unwrap(); + let mas_cb = |_: &mut CooMatrix| (); + let mut mm = CooMatrix::new(2, 2, 4, Sym::YesLower).unwrap(); + (mas_cb)(&mut mm); + system.set_mass(None, Sym::YesLower, mas_cb).unwrap(); + assert_eq!(system.get_mass_nnz(), 3); + system.set_mass(None, Sym::No, mas_cb).unwrap(); + assert_eq!(system.get_mass_nnz(), 4); + } } diff --git a/russell_ode/zscripts/memcheck.bash b/russell_ode/zscripts/memcheck.bash index 5f773aac..0667c4b3 100755 --- a/russell_ode/zscripts/memcheck.bash +++ b/russell_ode/zscripts/memcheck.bash @@ -4,7 +4,7 @@ INTEL_MKL=${1:-""} # 0 or 1 to use intel_mkl FEAT="" if [ "${INTEL_MKL}" = "1" ]; then - FEAT="--features intel_mkl" + FEAT="--features intel_mkl,local_suitesparse" fi cargo build $FEAT From c546e326f73f5341e316f1b854f68d7871b577c5 Mon Sep 17 00:00:00 2001 From: Dorival Pedroso Date: Wed, 22 May 2024 20:17:52 +1000 Subject: [PATCH 2/2] [Important] Impl a move of System into OdeSolver instead of taking a reference. Make System cloneable --- README.md | 2 +- russell_ode/README.md | 4 +- russell_ode/examples/amplifier1t_radau5.rs | 2 +- russell_ode/examples/arenstorf_dopri8.rs | 2 +- .../examples/brusselator_ode_dopri8.rs | 2 +- .../examples/brusselator_ode_fix_step.rs | 2 +- .../examples/brusselator_ode_var_step.rs | 2 +- .../brusselator_pde_2nd_comparison.rs | 2 +- .../examples/brusselator_pde_radau5.rs | 2 +- .../examples/brusselator_pde_radau5_2nd.rs | 2 +- russell_ode/examples/hairer_wanner_eq1.rs | 7 +-- .../pde_1d_heat_spectral_collocation.rs | 2 +- russell_ode/examples/robertson.rs | 4 +- .../examples/simple_ode_single_equation.rs | 2 +- .../examples/simple_system_with_mass.rs | 2 +- russell_ode/examples/van_der_pol_dopri5.rs | 2 +- russell_ode/examples/van_der_pol_radau5.rs | 2 +- russell_ode/src/bin/amplifier1t.rs | 2 +- russell_ode/src/bin/brusselator_pde.rs | 8 +-- russell_ode/src/euler_backward.rs | 17 ++++--- russell_ode/src/euler_forward.rs | 8 +-- russell_ode/src/explicit_runge_kutta.rs | 22 ++++----- russell_ode/src/lib.rs | 2 +- russell_ode/src/ode_solver.rs | 49 ++++++++++--------- russell_ode/src/output.rs | 5 -- russell_ode/src/radau5.rs | 23 ++++----- russell_ode/src/samples.rs | 28 +++++------ russell_ode/src/system.rs | 41 +++++++++++----- russell_ode/tests/test_bweuler.rs | 6 +-- russell_ode/tests/test_dopri5_arenstorf.rs | 2 +- .../tests/test_dopri5_arenstorf_debug.rs | 2 +- .../tests/test_dopri5_hairer_wanner_eq1.rs | 2 +- .../tests/test_dopri5_van_der_pol_debug.rs | 2 +- russell_ode/tests/test_dopri8_van_der_pol.rs | 2 +- .../tests/test_dopri8_van_der_pol_debug.rs | 2 +- russell_ode/tests/test_fweuler.rs | 2 +- russell_ode/tests/test_mdeuler.rs | 2 +- russell_ode/tests/test_radau5_amplifier1t.rs | 2 +- .../tests/test_radau5_brusselator_pde.rs | 2 +- .../tests/test_radau5_hairer_wanner_eq1.rs | 2 +- .../test_radau5_hairer_wanner_eq1_debug.rs | 2 +- russell_ode/tests/test_radau5_robertson.rs | 2 +- .../tests/test_radau5_robertson_debug.rs | 2 +- .../tests/test_radau5_robertson_small_h.rs | 2 +- russell_ode/tests/test_radau5_van_der_pol.rs | 2 +- .../tests/test_radau5_van_der_pol_debug.rs | 2 +- 46 files changed, 154 insertions(+), 134 deletions(-) diff --git a/README.md b/README.md index ab1236a4..64859ed3 100644 --- a/README.md +++ b/README.md @@ -597,7 +597,7 @@ fn main() -> Result<(), StrError> { // solver let params = Params::new(Method::DoPri8); - let mut solver = OdeSolver::new(params, &system)?; + let mut solver = OdeSolver::new(params, system)?; // enable dense output let h_out = 0.01; diff --git a/russell_ode/README.md b/russell_ode/README.md index 4642dd19..42ea2592 100644 --- a/russell_ode/README.md +++ b/russell_ode/README.md @@ -134,7 +134,7 @@ fn main() -> Result<(), StrError> { // solver let params = Params::new(Method::DoPri8); - let mut solver = OdeSolver::new(params, &system)?; + let mut solver = OdeSolver::new(params, system)?; // initial values let x = 0.0; @@ -269,7 +269,7 @@ fn main() -> Result<(), StrError> { // solver let params = Params::new(Method::Radau5); - let mut solver = OdeSolver::new(params, &system)?; + let mut solver = OdeSolver::new(params, system)?; // initial values let x = 0.0; diff --git a/russell_ode/examples/amplifier1t_radau5.rs b/russell_ode/examples/amplifier1t_radau5.rs index f1697235..1ecea77e 100644 --- a/russell_ode/examples/amplifier1t_radau5.rs +++ b/russell_ode/examples/amplifier1t_radau5.rs @@ -21,7 +21,7 @@ fn main() -> Result<(), StrError> { // solver let params = Params::new(Method::Radau5); - let mut solver = OdeSolver::new(params, &system)?; + let mut solver = OdeSolver::new(params, system)?; // enable dense output let h_out = 0.0001; diff --git a/russell_ode/examples/arenstorf_dopri8.rs b/russell_ode/examples/arenstorf_dopri8.rs index c5261b8d..b275576b 100644 --- a/russell_ode/examples/arenstorf_dopri8.rs +++ b/russell_ode/examples/arenstorf_dopri8.rs @@ -23,7 +23,7 @@ fn main() -> Result<(), StrError> { let params = Params::new(Method::DoPri8); // allocate the solver - let mut solver = OdeSolver::new(params, &system)?; + let mut solver = OdeSolver::new(params, system)?; // enable dense output let h_out = 0.01; diff --git a/russell_ode/examples/brusselator_ode_dopri8.rs b/russell_ode/examples/brusselator_ode_dopri8.rs index 41d1ab24..83062b2d 100644 --- a/russell_ode/examples/brusselator_ode_dopri8.rs +++ b/russell_ode/examples/brusselator_ode_dopri8.rs @@ -26,7 +26,7 @@ fn main() -> Result<(), StrError> { let params = Params::new(Method::DoPri8); // allocate the solver - let mut solver = OdeSolver::new(params, &system)?; + let mut solver = OdeSolver::new(params, system)?; // enable dense output let h_out = 0.01; diff --git a/russell_ode/examples/brusselator_ode_fix_step.rs b/russell_ode/examples/brusselator_ode_fix_step.rs index 5c15dc38..225fc4c5 100644 --- a/russell_ode/examples/brusselator_ode_fix_step.rs +++ b/russell_ode/examples/brusselator_ode_fix_step.rs @@ -44,7 +44,7 @@ fn main() -> Result<(), StrError> { for method in Method::erk_methods() { // allocate the solver let params = Params::new(method); - let mut solver = OdeSolver::new(params, &system)?; + let mut solver = OdeSolver::new(params, system.clone())?; // arrays holding the results let mut n_f_eval = vec![0.0; hh.len()]; diff --git a/russell_ode/examples/brusselator_ode_var_step.rs b/russell_ode/examples/brusselator_ode_var_step.rs index 65886e32..6618f2f5 100644 --- a/russell_ode/examples/brusselator_ode_var_step.rs +++ b/russell_ode/examples/brusselator_ode_var_step.rs @@ -44,7 +44,7 @@ fn main() -> Result<(), StrError> { for method in &[Method::Radau5, Method::Merson4, Method::DoPri5, Method::DoPri8] { // allocate the solver let params = Params::new(*method); - let mut solver = OdeSolver::new(params, &system)?; + let mut solver = OdeSolver::new(params, system.clone())?; // arrays holding the results let mut n_f_eval = vec![0.0; tols.len()]; diff --git a/russell_ode/examples/brusselator_pde_2nd_comparison.rs b/russell_ode/examples/brusselator_pde_2nd_comparison.rs index 77f038d1..f221a375 100644 --- a/russell_ode/examples/brusselator_pde_2nd_comparison.rs +++ b/russell_ode/examples/brusselator_pde_2nd_comparison.rs @@ -22,7 +22,7 @@ fn main() { params.set_tolerances(1e-4, 1e-4, None).unwrap(); // allocate the solver - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); // solve the ODE system let mut yy = yy0.clone(); diff --git a/russell_ode/examples/brusselator_pde_radau5.rs b/russell_ode/examples/brusselator_pde_radau5.rs index bc19d521..ac87064a 100644 --- a/russell_ode/examples/brusselator_pde_radau5.rs +++ b/russell_ode/examples/brusselator_pde_radau5.rs @@ -26,7 +26,7 @@ fn main() -> Result<(), StrError> { params.set_tolerances(1e-4, 1e-4, None)?; // allocate the solver - let mut solver = OdeSolver::new(params, &system)?; + let mut solver = OdeSolver::new(params, system)?; // output let h_out = 0.5; diff --git a/russell_ode/examples/brusselator_pde_radau5_2nd.rs b/russell_ode/examples/brusselator_pde_radau5_2nd.rs index a311b9f9..87cdbc34 100644 --- a/russell_ode/examples/brusselator_pde_radau5_2nd.rs +++ b/russell_ode/examples/brusselator_pde_radau5_2nd.rs @@ -27,7 +27,7 @@ fn main() -> Result<(), StrError> { params.set_tolerances(1e-4, 1e-4, None)?; // allocate the solver - let mut solver = OdeSolver::new(params, &system)?; + let mut solver = OdeSolver::new(params, system)?; // output let h_out = 1.0; diff --git a/russell_ode/examples/hairer_wanner_eq1.rs b/russell_ode/examples/hairer_wanner_eq1.rs index b4f4575f..b768047f 100644 --- a/russell_ode/examples/hairer_wanner_eq1.rs +++ b/russell_ode/examples/hairer_wanner_eq1.rs @@ -21,13 +21,14 @@ use russell_ode::prelude::*; fn main() -> Result<(), StrError> { // get the ODE system let (system, x0, y0, mut args, y_fn_x) = Samples::hairer_wanner_eq1(); + let ndim = system.get_ndim(); // final x let x1 = 1.5; // solvers - let mut bweuler = OdeSolver::new(Params::new(Method::BwEuler), &system)?; - let mut fweuler = OdeSolver::new(Params::new(Method::FwEuler), &system)?; + let mut bweuler = OdeSolver::new(Params::new(Method::BwEuler), system.clone())?; + let mut fweuler = OdeSolver::new(Params::new(Method::FwEuler), system)?; // solve the problem with BwEuler and h = 0.5 bweuler.enable_output().set_step_recording(&[0]); @@ -51,7 +52,7 @@ fn main() -> Result<(), StrError> { fweuler.solve(&mut y, x0, x1, Some(h), &mut args)?; // analytical solution - let mut y_aux = Vector::new(system.get_ndim()); + let mut y_aux = Vector::new(ndim); let mut x_ana1 = Vector::linspace(0.0, 0.2, 20)?; // small stepsizes let mut x_ana2 = Vector::linspace(0.201, x1, 20)?; // larger stepsizes x_ana1.as_mut_data().append(x_ana2.as_mut_data()); // merge the two vectors diff --git a/russell_ode/examples/pde_1d_heat_spectral_collocation.rs b/russell_ode/examples/pde_1d_heat_spectral_collocation.rs index c003905d..1f514b71 100644 --- a/russell_ode/examples/pde_1d_heat_spectral_collocation.rs +++ b/russell_ode/examples/pde_1d_heat_spectral_collocation.rs @@ -85,7 +85,7 @@ fn run( // ODE solver let mut params = Params::new(Method::DoPri8); params.set_tolerances(1e-10, 1e-10, None)?; - let mut ode = OdeSolver::new(params, &system)?; + let mut ode = OdeSolver::new(params, system)?; // interpolant let mut par = InterpParams::new(); diff --git a/russell_ode/examples/robertson.rs b/russell_ode/examples/robertson.rs index a1c51f99..a9afda50 100644 --- a/russell_ode/examples/robertson.rs +++ b/russell_ode/examples/robertson.rs @@ -43,8 +43,8 @@ fn main() -> Result<(), StrError> { params3.erk.lund_beta = 0.0; // solvers - let mut radau5 = OdeSolver::new(params1, &system)?; - let mut dopri5 = OdeSolver::new(params2, &system)?; + let mut radau5 = OdeSolver::new(params1, system.clone())?; + let mut dopri5 = OdeSolver::new(params2, system)?; // selected component for output let sel = 1; diff --git a/russell_ode/examples/simple_ode_single_equation.rs b/russell_ode/examples/simple_ode_single_equation.rs index 5cd5a9fb..03f53d32 100644 --- a/russell_ode/examples/simple_ode_single_equation.rs +++ b/russell_ode/examples/simple_ode_single_equation.rs @@ -11,7 +11,7 @@ fn main() -> Result<(), StrError> { // solver let params = Params::new(Method::DoPri8); - let mut solver = OdeSolver::new(params, &system)?; + let mut solver = OdeSolver::new(params, system)?; // initial values let x = 0.0; diff --git a/russell_ode/examples/simple_system_with_mass.rs b/russell_ode/examples/simple_system_with_mass.rs index ee2cca3d..7a7a05e0 100644 --- a/russell_ode/examples/simple_system_with_mass.rs +++ b/russell_ode/examples/simple_system_with_mass.rs @@ -40,7 +40,7 @@ fn main() -> Result<(), StrError> { // solver let params = Params::new(Method::Radau5); - let mut solver = OdeSolver::new(params, &system)?; + let mut solver = OdeSolver::new(params, system)?; // initial values let x = 0.0; diff --git a/russell_ode/examples/van_der_pol_dopri5.rs b/russell_ode/examples/van_der_pol_dopri5.rs index 40913546..28368430 100644 --- a/russell_ode/examples/van_der_pol_dopri5.rs +++ b/russell_ode/examples/van_der_pol_dopri5.rs @@ -30,7 +30,7 @@ fn main() -> Result<(), StrError> { params.set_tolerances(1e-5, 1e-5, None)?; // allocate the solver - let mut solver = OdeSolver::new(params, &system)?; + let mut solver = OdeSolver::new(params, system)?; // enable step and dense output let h_out = 0.01; diff --git a/russell_ode/examples/van_der_pol_radau5.rs b/russell_ode/examples/van_der_pol_radau5.rs index 9b07130c..e8660bdd 100644 --- a/russell_ode/examples/van_der_pol_radau5.rs +++ b/russell_ode/examples/van_der_pol_radau5.rs @@ -27,7 +27,7 @@ fn main() -> Result<(), StrError> { params.set_tolerances(1e-4, 1e-4, None)?; // allocate the solver - let mut solver = OdeSolver::new(params, &system)?; + let mut solver = OdeSolver::new(params, system)?; // enable step output let selected_y_components = &[0, 1]; diff --git a/russell_ode/src/bin/amplifier1t.rs b/russell_ode/src/bin/amplifier1t.rs index fd4f0526..6b8943a9 100644 --- a/russell_ode/src/bin/amplifier1t.rs +++ b/russell_ode/src/bin/amplifier1t.rs @@ -14,7 +14,7 @@ fn main() -> Result<(), StrError> { params.set_tolerances(1e-4, 1e-4, None)?; // allocate the solver - let mut solver = OdeSolver::new(params, &system)?; + let mut solver = OdeSolver::new(params, system)?; // enable dense output let h_out = 0.001; diff --git a/russell_ode/src/bin/brusselator_pde.rs b/russell_ode/src/bin/brusselator_pde.rs index 5f174410..deda93d9 100644 --- a/russell_ode/src/bin/brusselator_pde.rs +++ b/russell_ode/src/bin/brusselator_pde.rs @@ -99,7 +99,9 @@ fn main() -> Result<(), StrError> { } // solve the ODE system - let mut solver = OdeSolver::new(params, &system)?; + let ndim = system.get_ndim(); + let jac_nnz = system.get_jac_nnz(); + let mut solver = OdeSolver::new(params, system)?; solver.solve(&mut yy0, t0, t1, None, &mut args)?; // print stat @@ -108,8 +110,8 @@ fn main() -> Result<(), StrError> { println!("Number of points along x and y = {}", opt.npoint); println!("Tolerance (abs_tol = rel_tol) = {}", format_scientific(tol, 8, 2)); println!("Concurrent real and complex sys = {}", !opt.serial); - println!("Problem dimension (ndim) = {}", system.get_ndim()); - println!("Number of non-zeros (jac_nnz) = {}", system.get_jac_nnz()); + println!("Problem dimension (ndim) = {}", ndim); + println!("Number of non-zeros (jac_nnz) = {}", jac_nnz); println!("Number of BLAS threads = {}", get_num_threads()); println!("Linear solver = {:?}", params.newton.genie); println!("{}", stat); diff --git a/russell_ode/src/euler_backward.rs b/russell_ode/src/euler_backward.rs index b7067ac3..9923aca1 100644 --- a/russell_ode/src/euler_backward.rs +++ b/russell_ode/src/euler_backward.rs @@ -9,7 +9,7 @@ pub(crate) struct EulerBackward<'a, A> { params: Params, /// ODE system - system: &'a System<'a, A>, + system: System<'a, A>, /// Vector holding the function evaluation /// @@ -34,7 +34,7 @@ pub(crate) struct EulerBackward<'a, A> { impl<'a, A> EulerBackward<'a, A> { /// Allocates a new instance - pub fn new(params: Params, system: &'a System<'a, A>) -> Self { + pub fn new(params: Params, system: System<'a, A>) -> Self { let ndim = system.ndim; let jac_nnz = if params.newton.use_numerical_jacobian { ndim * ndim @@ -42,6 +42,7 @@ impl<'a, A> EulerBackward<'a, A> { system.jac_nnz }; let nnz = jac_nnz + ndim; // +ndim corresponds to the diagonal I matrix + let sym = system.symmetric; EulerBackward { params, system, @@ -49,7 +50,7 @@ impl<'a, A> EulerBackward<'a, A> { w: Vector::new(ndim), r: Vector::new(ndim), dy: Vector::new(ndim), - kk: SparseMatrix::new_coo(ndim, ndim, nnz, system.symmetric).unwrap(), + kk: SparseMatrix::new_coo(ndim, ndim, nnz, sym).unwrap(), solver: LinSolver::new(params.newton.genie).unwrap(), } } @@ -107,7 +108,7 @@ impl<'a, A> OdeSolverTrait for EulerBackward<'a, A> { work.stats.n_function += ndim; let w1 = &mut self.k; // workspace let w2 = &mut self.dy; // workspace - numerical_jacobian(kk, h, x_new, y_new, w1, w2, args, &self.system.function)?; + numerical_jacobian(kk, h, x_new, y_new, w1, w2, args, self.system.function.as_ref())?; } else { (self.system.jacobian.as_ref().unwrap())(kk, h, x_new, y_new, args)?; } @@ -260,7 +261,7 @@ mod tests { // allocate structs let params = Params::new(Method::BwEuler); - let mut solver = EulerBackward::new(params, &system); + let mut solver = EulerBackward::new(params, system); let mut work = Workspace::new(Method::BwEuler); // check dense output availability @@ -315,7 +316,7 @@ mod tests { // allocate structs let mut params = Params::new(Method::BwEuler); params.newton.use_numerical_jacobian = true; - let mut solver = EulerBackward::new(params, &system); + let mut solver = EulerBackward::new(params, system); let mut work = Workspace::new(Method::BwEuler); // numerical approximation @@ -359,7 +360,7 @@ mod tests { let mut params = Params::new(Method::BwEuler); let (system, x0, y0, mut args, _) = Samples::kreyszig_ex4_page920(); params.newton.n_iteration_max = 0; - let mut solver = EulerBackward::new(params, &system); + let mut solver = EulerBackward::new(params, system); let mut work = Workspace::new(Method::BwEuler); assert_eq!( solver.step(&mut work, x0, &y0, 0.1, &mut args).err(), @@ -391,7 +392,7 @@ mod tests { }) .unwrap(); let params = Params::new(Method::BwEuler); - let mut solver = EulerBackward::new(params, &system); + let mut solver = EulerBackward::new(params, system); let mut work = Workspace::new(Method::BwEuler); let x = 0.0; let y = Vector::from(&[0.0]); diff --git a/russell_ode/src/euler_forward.rs b/russell_ode/src/euler_forward.rs index dc8ba00d..67b645a2 100644 --- a/russell_ode/src/euler_forward.rs +++ b/russell_ode/src/euler_forward.rs @@ -8,7 +8,7 @@ use russell_lab::{vec_add, vec_copy, Vector}; /// and should not be used in production codes. pub(crate) struct EulerForward<'a, A> { /// ODE system - system: &'a System<'a, A>, + system: System<'a, A>, /// Vector holding the function evaluation /// @@ -21,7 +21,7 @@ pub(crate) struct EulerForward<'a, A> { impl<'a, A> EulerForward<'a, A> { /// Allocates a new instance - pub fn new(system: &'a System<'a, A>) -> Self { + pub fn new(system: System<'a, A>) -> Self { let ndim = system.ndim; EulerForward { system, @@ -86,7 +86,7 @@ mod tests { let ndim = system.ndim; // allocate structs - let mut solver = EulerForward::new(&system); + let mut solver = EulerForward::new(system); let mut work = Workspace::new(Method::FwEuler); // check dense output availability @@ -163,7 +163,7 @@ mod tests { f[0] = 1.0; Err("stop") }); - let mut solver = EulerForward::new(&system); + let mut solver = EulerForward::new(system); let mut work = Workspace::new(Method::FwEuler); let x = 0.0; let y = Vector::from(&[0.0]); diff --git a/russell_ode/src/explicit_runge_kutta.rs b/russell_ode/src/explicit_runge_kutta.rs index 061e9885..ca125e91 100644 --- a/russell_ode/src/explicit_runge_kutta.rs +++ b/russell_ode/src/explicit_runge_kutta.rs @@ -25,7 +25,7 @@ pub(crate) struct ExplicitRungeKutta<'a, A> { params: Params, /// ODE system - system: &'a System<'a, A>, + system: System<'a, A>, /// Information such as implicit, embedded, etc. info: Information, @@ -77,7 +77,7 @@ pub(crate) struct ExplicitRungeKutta<'a, A> { impl<'a, A> ExplicitRungeKutta<'a, A> { /// Allocates a new instance - pub fn new(params: Params, system: &'a System<'a, A>) -> Result { + pub fn new(params: Params, system: System<'a, A>) -> Result { // Runge-Kutta coefficients #[rustfmt::skip] let (aa, bb, cc) = match params.method { @@ -377,7 +377,7 @@ mod tests { println!("\n... {:?} ...", method); let params = Params::new(method); let system = System::new(1, |_, _, _, _args: &mut Args| Ok(())); - let erk = ExplicitRungeKutta::new(params, &system).unwrap(); + let erk = ExplicitRungeKutta::new(params, system).unwrap(); let nstage = erk.nstage; assert_eq!(erk.aa.dims(), (nstage, nstage)); assert_eq!(erk.bb.dim(), nstage); @@ -502,7 +502,7 @@ mod tests { // allocate structs let params = Params::new(Method::MdEuler); // aka the Improved Euler in Kreyszig's book - let mut solver = ExplicitRungeKutta::new(params, &system).unwrap(); + let mut solver = ExplicitRungeKutta::new(params, system).unwrap(); let mut work = Workspace::new(Method::FwEuler); // check dense output availability @@ -585,7 +585,7 @@ mod tests { // allocate structs let params = Params::new(Method::Rk4); // aka the Classical RK in Kreyszig's book - let mut solver = ExplicitRungeKutta::new(params, &system).unwrap(); + let mut solver = ExplicitRungeKutta::new(params, system).unwrap(); let mut work = Workspace::new(Method::FwEuler); // check dense output availability @@ -686,7 +686,7 @@ mod tests { // allocate solver let params = Params::new(Method::Fehlberg4); - let mut solver = ExplicitRungeKutta::new(params, &system).unwrap(); + let mut solver = ExplicitRungeKutta::new(params, system).unwrap(); let mut work = Workspace::new(Method::FwEuler); // perform one step (compute k) @@ -733,20 +733,20 @@ mod tests { }); assert_eq!( - ExplicitRungeKutta::new(Params::new(Method::Radau5), &system).err(), + ExplicitRungeKutta::new(Params::new(Method::Radau5), system.clone()).err(), Some("cannot use Radau5 with ExplicitRungeKutta") ); assert_eq!( - ExplicitRungeKutta::new(Params::new(Method::BwEuler), &system).err(), + ExplicitRungeKutta::new(Params::new(Method::BwEuler), system.clone()).err(), Some("cannot use BwEuler with ExplicitRungeKutta") ); assert_eq!( - ExplicitRungeKutta::new(Params::new(Method::FwEuler), &system).err(), + ExplicitRungeKutta::new(Params::new(Method::FwEuler), system.clone()).err(), Some("cannot use FwEuler with ExplicitRungeKutta") ); let params = Params::new(Method::DoPri8); - let mut solver = ExplicitRungeKutta::new(params, &system).unwrap(); + let mut solver = ExplicitRungeKutta::new(params, system).unwrap(); let mut work = Workspace::new(Method::DoPri8); let mut x = 0.0; let mut y = Vector::from(&[0.0]); @@ -782,7 +782,7 @@ mod tests { for method in Method::erk_methods() { let params = Params::new(method); let mut work = Workspace::new(method); - let mut solver = ExplicitRungeKutta::new(params, &system).unwrap(); + let mut solver = ExplicitRungeKutta::new(params, system.clone()).unwrap(); solver.step(&mut work, x, &y, h, &mut args).unwrap(); work.stats.n_accepted += 1; // important (must precede accept) solver.accept(&mut work, &mut x, &mut y, h, &mut args).unwrap(); diff --git a/russell_ode/src/lib.rs b/russell_ode/src/lib.rs index f3e1047b..b3749821 100644 --- a/russell_ode/src/lib.rs +++ b/russell_ode/src/lib.rs @@ -172,7 +172,7 @@ //! //! // solver //! let params = Params::new(Method::Radau5); -//! let mut solver = OdeSolver::new(params, &system)?; +//! let mut solver = OdeSolver::new(params, system)?; //! //! // initial values //! let x = 0.0; diff --git a/russell_ode/src/ode_solver.rs b/russell_ode/src/ode_solver.rs index 963a93ad..20221b3e 100644 --- a/russell_ode/src/ode_solver.rs +++ b/russell_ode/src/ode_solver.rs @@ -83,7 +83,7 @@ use russell_lab::{vec_all_finite, Vector}; /// /// // solver /// let params = Params::new(Method::Radau5); -/// let mut solver = OdeSolver::new(params, &system)?; +/// let mut solver = OdeSolver::new(params, system)?; /// /// // initial values /// let x = 0.0; @@ -135,14 +135,14 @@ impl<'a, A> OdeSolver<'a, A> { /// /// * `A` -- generic argument to assist in the f(x,y) and Jacobian functions. /// It may be simply [NoArgs] indicating that no arguments are needed. - pub fn new(params: Params, system: &'a System<'a, A>) -> Result + pub fn new(params: Params, system: System<'a, A>) -> Result where A: 'a, { + params.validate()?; if system.calc_mass.is_some() && params.method != Method::Radau5 { return Err("the method must be Radau5 for systems with a mass matrix"); } - params.validate()?; let ndim = system.ndim; let actual: Box> = if params.method == Method::Radau5 { Box::new(Radau5::new(params, system)) @@ -478,13 +478,13 @@ mod tests { let (system, _, _, _, _) = Samples::simple_system_with_mass_matrix(false, Genie::Umfpack); let mut params = Params::new(Method::MdEuler); assert_eq!( - OdeSolver::new(params, &system).err(), + OdeSolver::new(params, system).err(), Some("the method must be Radau5 for systems with a mass matrix") ); let (system, _, _, _, _) = Samples::simple_equation_constant(); params.step.m_max = 0.0; // wrong assert_eq!( - OdeSolver::new(params, &system).err(), + OdeSolver::new(params, system).err(), Some("parameter must satisfy: 0.001 ≤ m_min < 0.5 and m_min < m_max") ); } @@ -492,14 +492,15 @@ mod tests { #[test] fn solve_captures_errors() { let (system, _, _, mut args, _) = Samples::simple_equation_constant(); + let ndim = system.ndim; let params = Params::new(Method::FwEuler); - let mut solver = OdeSolver::new(params, &system).unwrap(); - let mut y0 = Vector::new(system.ndim + 1); // wrong dim + let mut solver = OdeSolver::new(params, system).unwrap(); + let mut y0 = Vector::new(ndim + 1); // wrong dim assert_eq!( solver.solve(&mut y0, 0.0, 1.0, None, &mut args).err(), Some("y0.dim() must be equal to ndim") ); - let mut y0 = Vector::new(system.ndim); + let mut y0 = Vector::new(ndim); assert_eq!( solver.solve(&mut y0, 0.0, 0.0, None, &mut args).err(), Some("x1 must be greater than x0") @@ -517,13 +518,13 @@ mod tests { // it becomes stiff and yield infinite results let (system, _, mut y0, mut args, _) = Samples::brusselator_ode(); let params = Params::new(Method::FwEuler); - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system.clone()).unwrap(); assert_eq!( solver.solve(&mut y0, 0.0, 9.0, Some(1.0), &mut args).err(), Some("an element of the vector is either infinite or NaN") ); let params = Params::new(Method::MdEuler); - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); assert_eq!( solver.solve(&mut y0, 0.0, 1.0, None, &mut args).err(), Some("an element of the vector is either infinite or NaN") @@ -535,7 +536,7 @@ mod tests { let (system, _, mut y0, mut args, _) = Samples::simple_equation_constant(); let mut params = Params::new(Method::MdEuler); params.step.n_step_max = 1; // will make the solver to fail (too few steps) - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); assert_eq!( solver.solve(&mut y0, 0.0, 1.0, None, &mut args).err(), Some("variable stepping did not converge") @@ -546,7 +547,7 @@ mod tests { fn out_initialize_errors_are_captured() { let (system, _, mut y0, mut args, _) = Samples::simple_equation_constant(); let params = Params::new(Method::DoPri5); - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); solver .enable_output() .set_dense_x_out(&[0.0, 1.0]) @@ -569,7 +570,7 @@ mod tests { let (system, x0, y0, mut args, _) = Samples::simple_equation_constant(); let x1 = 1.0; let params = Params::new(Method::FwEuler); - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); let mut y = y0.clone(); solver.solve(&mut y, x0, x1, None, &mut args).unwrap(); vec_approx_eq(&y, &[1.0], 1e-15); @@ -581,7 +582,7 @@ mod tests { let (system, _, mut y0, mut args, _) = Samples::simple_equation_constant(); let mut params = Params::new(Method::DoPri5); params.step.h_ini = 20.0; // will be truncated to 1 yielding a single step - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); solver.solve(&mut y0, 0.0, 1.0, None, &mut args).unwrap(); assert_eq!(solver.work.stats.n_accepted, 1); vec_approx_eq(&y0, &[1.0], 1e-15); @@ -592,7 +593,7 @@ mod tests { let (system, _, mut y0, mut args, _) = Samples::simple_equation_constant(); let mut params = Params::new(Method::MdEuler); params.step.h_ini = 0.1; - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); solver.solve(&mut y0, 0.0, 0.3, None, &mut args).unwrap(); vec_approx_eq(&y0, &[0.3], 1e-15); } @@ -603,11 +604,11 @@ mod tests { let mut params = Params::new(Method::MdEuler); params.step.n_step_max = 0; assert_eq!( - OdeSolver::new(params, &system).err(), + OdeSolver::new(params, system.clone()).err(), Some("parameter must satisfy: n_step_max ≥ 1") ); params.step.n_step_max = 1000; - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); assert_eq!(solver.params.step.n_step_max, 1000); params.step.n_step_max = 2; // this will not change the solver until update_params is called assert_eq!(solver.params.step.n_step_max, 1000); @@ -631,7 +632,7 @@ mod tests { let (system, x0, mut y0, mut args, _) = Samples::simple_equation_constant(); let x1 = 1.0; let params = Params::new(Method::FwEuler); - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); solver.enable_output().set_dense_recording(&[0]); assert_eq!( solver.solve(&mut y0, x0, x1, None, &mut args).err(), @@ -644,7 +645,7 @@ mod tests { // system and solver let (system, _, y0, mut args, y_fn_x) = Samples::simple_equation_constant(); let params = Params::new(Method::DoPri5); - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); // output let path_key = "/tmp/russell_ode/test_solve_step_output_works"; @@ -708,7 +709,7 @@ mod tests { // system and solver let (system, _, y0, mut args, _) = Samples::simple_equation_constant(); let params = Params::new(Method::FwEuler); - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); // output solver.enable_output().set_step_recording(&[0]); @@ -743,7 +744,7 @@ mod tests { // system and solver let (system, _, _, mut args, _) = Samples::simple_equation_constant(); let params = Params::new(Method::DoPri5); - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); // output solver @@ -771,7 +772,7 @@ mod tests { // system and solver let (system, _, y0, mut args, y_fn_x) = Samples::simple_equation_constant(); let params = Params::new(Method::DoPri5); - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); // output const H_OUT: f64 = 0.1; @@ -874,7 +875,7 @@ mod tests { // system and solver let (system, _, _, mut args, _) = Samples::simple_equation_constant(); let params = Params::new(Method::DoPri5); - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); // interior output stations and selected y component let interior_x_out = &[-0.5, 0.0, 0.5]; @@ -936,7 +937,7 @@ mod tests { // parameters and solver let mut params = Params::new(Method::DoPri8); params.step.h_ini = 0.2; - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); // output solver.enable_output().set_dense_h_out(0.1).unwrap().set_dense_callback( diff --git a/russell_ode/src/output.rs b/russell_ode/src/output.rs index 22a1eb14..50bd7a28 100644 --- a/russell_ode/src/output.rs +++ b/russell_ode/src/output.rs @@ -5,7 +5,6 @@ use serde::{Deserialize, Serialize}; use std::collections::HashMap; use std::fs::{self, File}; use std::io::BufReader; -use std::marker::PhantomData; use std::path::Path; /// Holds the data generated at an accepted step or during the dense output @@ -125,9 +124,6 @@ pub struct Output<'a, A> { /// Holds the y(x) function (e.g., to compute the correct/analytical solution) yx_function: Option>, - - /// Handles the generic argument - phantom: PhantomData A>, } impl OutData { @@ -210,7 +206,6 @@ impl<'a, A> Output<'a, A> { // auxiliary y_aux: Vector::new(EMPTY), yx_function: None, - phantom: PhantomData, } } diff --git a/russell_ode/src/radau5.rs b/russell_ode/src/radau5.rs index 8663ac97..d46d3980 100644 --- a/russell_ode/src/radau5.rs +++ b/russell_ode/src/radau5.rs @@ -29,7 +29,7 @@ pub(crate) struct Radau5<'a, A> { params: Params, /// ODE system - system: &'a System<'a, A>, + system: System<'a, A>, /// Holds the mass matrix mass: Option, @@ -120,7 +120,7 @@ pub(crate) struct Radau5<'a, A> { impl<'a, A> Radau5<'a, A> { /// Allocates a new instance - pub fn new(params: Params, system: &'a System<'a, A>) -> Self { + pub fn new(params: Params, system: System<'a, A>) -> Self { let ndim = system.ndim; let (mass, mass_nnz) = match system.calc_mass.as_ref() { Some(calc) => { @@ -140,14 +140,15 @@ impl<'a, A> Radau5<'a, A> { system.jac_nnz }; let nnz = mass_nnz + jac_nnz; + let sym = system.symmetric; let theta = params.radau5.theta_max; Radau5 { params, system, mass, - jj: SparseMatrix::new_coo(ndim, ndim, jac_nnz, system.symmetric).unwrap(), - kk_real: SparseMatrix::new_coo(ndim, ndim, nnz, system.symmetric).unwrap(), - kk_comp: ComplexSparseMatrix::new_coo(ndim, ndim, nnz, system.symmetric).unwrap(), + jj: SparseMatrix::new_coo(ndim, ndim, jac_nnz, sym).unwrap(), + kk_real: SparseMatrix::new_coo(ndim, ndim, nnz, sym).unwrap(), + kk_comp: ComplexSparseMatrix::new_coo(ndim, ndim, nnz, sym).unwrap(), solver_real: LinSolver::new(params.newton.genie).unwrap(), solver_comp: ComplexLinSolver::new(params.newton.genie).unwrap(), reuse_jacobian: false, @@ -208,7 +209,7 @@ impl<'a, A> Radau5<'a, A> { let w1 = &mut self.dw0; // workspace let w2 = &mut self.dw1; // workspace vec_copy(y_mut, y).unwrap(); - numerical_jacobian(jj, 1.0, x, y_mut, w1, w2, args, &self.system.function)?; + numerical_jacobian(jj, 1.0, x, y_mut, w1, w2, args, self.system.function.as_ref())?; } else { (self.system.jacobian.as_ref().unwrap())(jj, 1.0, x, y, args)?; } @@ -756,7 +757,7 @@ mod tests { // allocate structs let params = Params::new(Method::Radau5); - let mut solver = Radau5::new(params, &system); + let mut solver = Radau5::new(params, system); let mut work = Workspace::new(Method::Radau5); // message @@ -825,7 +826,7 @@ mod tests { // allocate structs let mut params = Params::new(Method::Radau5); params.newton.use_numerical_jacobian = true; - let mut solver = Radau5::new(params, &system); + let mut solver = Radau5::new(params, system); let mut work = Workspace::new(Method::Radau5); // message @@ -907,7 +908,7 @@ mod tests { }) .unwrap(); let params = Params::new(Method::Radau5); - let mut solver = Radau5::new(params, &system); + let mut solver = Radau5::new(params, system); let mut work = Workspace::new(Method::Radau5); let x = 0.0; let y = Vector::from(&[0.0]); @@ -940,7 +941,7 @@ mod tests { let mut params = Params::new(Method::Radau5); params.newton.genie = genie; params.newton.use_numerical_jacobian = numerical; - let mut solver = Radau5::new(params, &system); + let mut solver = Radau5::new(params, system); let mut work = Workspace::new(Method::Radau5); // message @@ -1004,7 +1005,7 @@ mod tests { let mut params = Params::new(Method::Radau5); params.newton.genie = genie; params.newton.use_numerical_jacobian = numerical; - let mut solver = Radau5::new(params, &system); + let mut solver = Radau5::new(params, system); let mut work = Workspace::new(Method::Radau5); // message diff --git a/russell_ode/src/samples.rs b/russell_ode/src/samples.rs index 9a8b1fd0..0c9edebd 100644 --- a/russell_ode/src/samples.rs +++ b/russell_ode/src/samples.rs @@ -1286,7 +1286,7 @@ mod tests { (jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); // compute the numerical Jacobian matrix - let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function).unwrap(); + let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function.as_ref()).unwrap(); // check the Jacobian matrix let ana = jj.as_dense(); @@ -1312,7 +1312,7 @@ mod tests { (jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); // compute the numerical Jacobian matrix - let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function).unwrap(); + let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function.as_ref()).unwrap(); // check the Jacobian matrix let ana = jj.as_dense(); @@ -1332,7 +1332,7 @@ mod tests { (jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); // compute the numerical Jacobian matrix - let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function).unwrap(); + let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function.as_ref()).unwrap(); // check the Jacobian matrix let ana = jj.as_dense(); @@ -1354,7 +1354,7 @@ mod tests { (jacobian)(&mut jj, jac_alpha, x0, &y0, &mut args).unwrap(); // compute the numerical Jacobian matrix - let num = num_jacobian(system.ndim, x0, &y0, jac_alpha, &mut args, system.function).unwrap(); + let num = num_jacobian(system.ndim, x0, &y0, jac_alpha, &mut args, system.function.as_ref()).unwrap(); // check the Jacobian matrix let ana = jj.as_dense(); @@ -1376,7 +1376,7 @@ mod tests { (jacobian)(&mut jj, jac_alpha, x0, &y0, &mut args).unwrap(); // compute the numerical Jacobian matrix - let num = num_jacobian(system.ndim, x0, &y0, jac_alpha, &mut args, system.function).unwrap(); + let num = num_jacobian(system.ndim, x0, &y0, jac_alpha, &mut args, system.function.as_ref()).unwrap(); // check the Jacobian matrix let ana = jj.as_dense(); @@ -1401,7 +1401,7 @@ mod tests { (jacobian)(&mut jj, jac_alpha, t, &y0, &mut args).unwrap(); // compute the numerical Jacobian matrix - let num = num_jacobian(system.ndim, t, &y0, jac_alpha, &mut args, &system.function).unwrap(); + let num = num_jacobian(system.ndim, t, &y0, jac_alpha, &mut args, system.function.as_ref()).unwrap(); // check the Jacobian matrix let ana = jj.as_dense(); @@ -1422,7 +1422,7 @@ mod tests { (jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); // compute the numerical Jacobian matrix - let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function).unwrap(); + let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function.as_ref()).unwrap(); // check the Jacobian matrix let ana = jj.as_dense(); @@ -1448,7 +1448,7 @@ mod tests { (jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); // compute the numerical Jacobian matrix - let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function).unwrap(); + let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function.as_ref()).unwrap(); // check the Jacobian matrix let ana = jj.as_dense(); @@ -1468,7 +1468,7 @@ mod tests { (jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); // compute the numerical Jacobian matrix - let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function).unwrap(); + let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function.as_ref()).unwrap(); // check the Jacobian matrix let ana = jj.as_dense(); @@ -1488,7 +1488,7 @@ mod tests { (jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); // compute the numerical Jacobian matrix - let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function).unwrap(); + let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function.as_ref()).unwrap(); // check the Jacobian matrix let ana = jj.as_dense(); @@ -1508,7 +1508,7 @@ mod tests { (jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); // compute the numerical Jacobian matrix - let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function).unwrap(); + let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function.as_ref()).unwrap(); // check the Jacobian matrix let ana = jj.as_dense(); @@ -1528,7 +1528,7 @@ mod tests { (jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); // compute the numerical Jacobian matrix - let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function).unwrap(); + let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function.as_ref()).unwrap(); // check the Jacobian matrix let ana = jj.as_dense(); @@ -1563,7 +1563,7 @@ mod tests { (jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); // compute the numerical Jacobian matrix - let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function).unwrap(); + let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function.as_ref()).unwrap(); // check the Jacobian matrix let ana = jj.as_dense(); @@ -1589,7 +1589,7 @@ mod tests { (jacobian)(&mut jj, alpha, x0, &y0, &mut args).unwrap(); // compute the numerical Jacobian matrix - let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function).unwrap(); + let num = num_jacobian(system.ndim, x0, &y0, alpha, &mut args, system.function.as_ref()).unwrap(); // check the Jacobian matrix let ana = jj.as_dense(); diff --git a/russell_ode/src/system.rs b/russell_ode/src/system.rs index 6fa93e59..6e5e5b2e 100644 --- a/russell_ode/src/system.rs +++ b/russell_ode/src/system.rs @@ -1,7 +1,7 @@ use crate::StrError; use russell_lab::Vector; use russell_sparse::{CooMatrix, Sym}; -use std::marker::PhantomData; +use std::sync::Arc; /// Indicates that the system functions do not require extra arguments pub type NoArgs = u8; @@ -66,13 +66,13 @@ pub struct System<'a, A> { pub(crate) ndim: usize, /// ODE system function - pub(crate) function: Box Result<(), StrError> + 'a>, + pub(crate) function: Arc Result<(), StrError> + 'a>, /// Jacobian function - pub(crate) jacobian: Option Result<(), StrError> + 'a>>, + pub(crate) jacobian: Option Result<(), StrError> + 'a>>, /// Calc mass matrix function - pub(crate) calc_mass: Option>, + pub(crate) calc_mass: Option>, /// Number of non-zeros in the Jacobian matrix pub(crate) jac_nnz: usize, @@ -88,9 +88,6 @@ pub struct System<'a, A> { /// Symmetric type of the Jacobian and mass matrices pub(crate) symmetric: Sym, - - /// Handle generic argument - phantom: PhantomData A>, } impl<'a, A> System<'a, A> { @@ -157,7 +154,7 @@ impl<'a, A> System<'a, A> { pub fn new(ndim: usize, function: impl Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError> + 'a) -> Self { System { ndim, - function: Box::new(function), + function: Arc::new(function), jacobian: None, calc_mass: None, jac_nnz: ndim * ndim, @@ -165,7 +162,21 @@ impl<'a, A> System<'a, A> { sym_jac: None, sym_mass: None, symmetric: Sym::No, - phantom: PhantomData, + } + } + + /// Returns a copy of this struct + pub fn clone(&self) -> Self { + System { + ndim: self.ndim, + function: self.function.clone(), + jacobian: self.jacobian.clone(), + calc_mass: self.calc_mass.clone(), + jac_nnz: self.jac_nnz, + mass_nnz: self.mass_nnz, + sym_jac: self.sym_jac, + sym_mass: self.sym_mass, + symmetric: self.symmetric, } } @@ -202,7 +213,7 @@ impl<'a, A> System<'a, A> { }; self.sym_jac = Some(symmetric); self.symmetric = symmetric; - self.jacobian = Some(Box::new(callback)); + self.jacobian = Some(Arc::new(callback)); Ok(()) } @@ -237,7 +248,7 @@ impl<'a, A> System<'a, A> { }; self.sym_mass = Some(symmetric); self.symmetric = symmetric; - self.calc_mass = Some(Box::new(callback)); + self.calc_mass = Some(Arc::new(callback)); Ok(()) } @@ -328,6 +339,14 @@ mod tests { println!("n_function_eval = {}", args.n_function_eval); assert_eq!(args.n_function_eval, 1); assert_eq!(args.more_data_goes_here, true); + // clone + let clone = system.clone(); + assert_eq!(clone.ndim, 2); + assert_eq!(clone.jac_nnz, 4); + assert_eq!(clone.mass_nnz, 0); + assert_eq!(clone.sym_jac, None); + assert_eq!(clone.sym_mass, None); + assert_eq!(clone.symmetric, Sym::No); } #[test] diff --git a/russell_ode/tests/test_bweuler.rs b/russell_ode/tests/test_bweuler.rs index 227d52d0..c9c1469b 100644 --- a/russell_ode/tests/test_bweuler.rs +++ b/russell_ode/tests/test_bweuler.rs @@ -12,7 +12,7 @@ fn test_bweuler_hairer_wanner_eq1() { // set configuration parameters let params = Params::new(Method::BwEuler); - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); // solve the ODE system let h_equal = Some(1.875 / 50.0); @@ -57,7 +57,7 @@ fn test_bweuler_hairer_wanner_eq1_num_jac() { params.newton.use_numerical_jacobian = true; // solve the ODE system - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); let h_equal = Some(1.875 / 50.0); solver.solve(&mut y0, x0, x1, h_equal, &mut args).unwrap(); @@ -100,7 +100,7 @@ fn test_bweuler_hairer_wanner_eq1_modified_newton() { params.bweuler.use_modified_newton = true; // solve the ODE system - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); let h_equal = Some(1.875 / 50.0); solver.solve(&mut y0, x0, x1, h_equal, &mut args).unwrap(); diff --git a/russell_ode/tests/test_dopri5_arenstorf.rs b/russell_ode/tests/test_dopri5_arenstorf.rs index c8fc14d7..29a224d2 100644 --- a/russell_ode/tests/test_dopri5_arenstorf.rs +++ b/russell_ode/tests/test_dopri5_arenstorf.rs @@ -12,7 +12,7 @@ fn test_dopri5_arenstorf() { params.set_tolerances(1e-7, 1e-7, None).unwrap(); // allocate the solver - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); // enable dense output with 1.0 spacing solver diff --git a/russell_ode/tests/test_dopri5_arenstorf_debug.rs b/russell_ode/tests/test_dopri5_arenstorf_debug.rs index b2ca3a71..a2f32ebb 100644 --- a/russell_ode/tests/test_dopri5_arenstorf_debug.rs +++ b/russell_ode/tests/test_dopri5_arenstorf_debug.rs @@ -13,7 +13,7 @@ fn test_dopri5_arenstorf_debug() { params.debug = true; // allocate the solver - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); // solve the ODE system let y = &mut y0; diff --git a/russell_ode/tests/test_dopri5_hairer_wanner_eq1.rs b/russell_ode/tests/test_dopri5_hairer_wanner_eq1.rs index f731c6c3..58fc6857 100644 --- a/russell_ode/tests/test_dopri5_hairer_wanner_eq1.rs +++ b/russell_ode/tests/test_dopri5_hairer_wanner_eq1.rs @@ -15,7 +15,7 @@ fn test_dopri5_hairer_wanner_eq1() { params.step.h_ini = 1e-4; // allocate the solver - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); // enable dense output solver diff --git a/russell_ode/tests/test_dopri5_van_der_pol_debug.rs b/russell_ode/tests/test_dopri5_van_der_pol_debug.rs index 29eef983..16b822a7 100644 --- a/russell_ode/tests/test_dopri5_van_der_pol_debug.rs +++ b/russell_ode/tests/test_dopri5_van_der_pol_debug.rs @@ -18,7 +18,7 @@ fn test_dopri5_van_der_pol_debug() { params.stiffness.save_results = true; // allocate the solver - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); // enable output (to save stiff stations) solver.enable_output(); diff --git a/russell_ode/tests/test_dopri8_van_der_pol.rs b/russell_ode/tests/test_dopri8_van_der_pol.rs index c78ddc86..3cf84525 100644 --- a/russell_ode/tests/test_dopri8_van_der_pol.rs +++ b/russell_ode/tests/test_dopri8_van_der_pol.rs @@ -13,7 +13,7 @@ fn test_dopri8_van_der_pol() { params.set_tolerances(1e-9, 1e-9, None).unwrap(); // allocate the solver - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); // enable dense output solver diff --git a/russell_ode/tests/test_dopri8_van_der_pol_debug.rs b/russell_ode/tests/test_dopri8_van_der_pol_debug.rs index 2b9b1be6..b98012d0 100644 --- a/russell_ode/tests/test_dopri8_van_der_pol_debug.rs +++ b/russell_ode/tests/test_dopri8_van_der_pol_debug.rs @@ -18,7 +18,7 @@ fn test_dopri8_van_der_pol_debug() { params.stiffness.save_results = true; // allocate the solver - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); // output (to save stiff stations) solver.enable_output(); diff --git a/russell_ode/tests/test_fweuler.rs b/russell_ode/tests/test_fweuler.rs index 8383ca8d..d97a788b 100644 --- a/russell_ode/tests/test_fweuler.rs +++ b/russell_ode/tests/test_fweuler.rs @@ -14,7 +14,7 @@ fn test_fweuler_hairer_wanner_eq1() { let params = Params::new(Method::FwEuler); // solve the ODE system - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); let h_equal = Some(1.875 / 50.0); solver.solve(&mut y0, x0, x1, h_equal, &mut args).unwrap(); diff --git a/russell_ode/tests/test_mdeuler.rs b/russell_ode/tests/test_mdeuler.rs index ec31ee61..71c55404 100644 --- a/russell_ode/tests/test_mdeuler.rs +++ b/russell_ode/tests/test_mdeuler.rs @@ -15,7 +15,7 @@ fn test_mdeuler_hairer_wanner_eq1() { params.step.h_ini = 1e-4; // solve the ODE system - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); solver.solve(&mut y0, x0, x1, None, &mut args).unwrap(); // get statistics diff --git a/russell_ode/tests/test_radau5_amplifier1t.rs b/russell_ode/tests/test_radau5_amplifier1t.rs index 98a7bf2e..81eff975 100644 --- a/russell_ode/tests/test_radau5_amplifier1t.rs +++ b/russell_ode/tests/test_radau5_amplifier1t.rs @@ -15,7 +15,7 @@ fn test_radau5_amplifier1t() { params.set_tolerances(1e-4, 1e-4, None).unwrap(); // allocate the solver - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); // enable output of accepted steps solver diff --git a/russell_ode/tests/test_radau5_brusselator_pde.rs b/russell_ode/tests/test_radau5_brusselator_pde.rs index af070951..cee9d5f4 100644 --- a/russell_ode/tests/test_radau5_brusselator_pde.rs +++ b/russell_ode/tests/test_radau5_brusselator_pde.rs @@ -21,7 +21,7 @@ fn test_radau5_brusselator_pde() { params.set_tolerances(1e-3, 1e-3, None).unwrap(); // allocate the solver - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); // solve the ODE system let yy = &mut yy0; diff --git a/russell_ode/tests/test_radau5_hairer_wanner_eq1.rs b/russell_ode/tests/test_radau5_hairer_wanner_eq1.rs index 5791c6cf..5b488bfa 100644 --- a/russell_ode/tests/test_radau5_hairer_wanner_eq1.rs +++ b/russell_ode/tests/test_radau5_hairer_wanner_eq1.rs @@ -15,7 +15,7 @@ fn test_radau5_hairer_wanner_eq1() { params.step.h_ini = 1e-4; // allocate the solver - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); // enable dense output solver diff --git a/russell_ode/tests/test_radau5_hairer_wanner_eq1_debug.rs b/russell_ode/tests/test_radau5_hairer_wanner_eq1_debug.rs index 027869ac..8bdf32d1 100644 --- a/russell_ode/tests/test_radau5_hairer_wanner_eq1_debug.rs +++ b/russell_ode/tests/test_radau5_hairer_wanner_eq1_debug.rs @@ -16,7 +16,7 @@ fn test_radau5_hairer_wanner_eq1_debug() { params.debug = true; // allocate the solver - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); // solve the ODE system solver.solve(&mut y0, x0, x1, None, &mut args).unwrap(); diff --git a/russell_ode/tests/test_radau5_robertson.rs b/russell_ode/tests/test_radau5_robertson.rs index 5adc179c..56e1788f 100644 --- a/russell_ode/tests/test_radau5_robertson.rs +++ b/russell_ode/tests/test_radau5_robertson.rs @@ -15,7 +15,7 @@ fn test_radau5_robertson() { params.set_tolerances(1e-8, 1e-2, None).unwrap(); // allocate the solver - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); // enable output of accepted steps solver.enable_output().set_step_recording(&[0, 1, 2]); diff --git a/russell_ode/tests/test_radau5_robertson_debug.rs b/russell_ode/tests/test_radau5_robertson_debug.rs index 21ff0fb5..73ad9ced 100644 --- a/russell_ode/tests/test_radau5_robertson_debug.rs +++ b/russell_ode/tests/test_radau5_robertson_debug.rs @@ -16,7 +16,7 @@ fn test_radau5_robertson_debug() { params.debug = true; // allocate the solver - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); // solve the ODE system solver.solve(&mut y0, x0, x1, None, &mut args).unwrap(); diff --git a/russell_ode/tests/test_radau5_robertson_small_h.rs b/russell_ode/tests/test_radau5_robertson_small_h.rs index f8882868..bc5f38ce 100644 --- a/russell_ode/tests/test_radau5_robertson_small_h.rs +++ b/russell_ode/tests/test_radau5_robertson_small_h.rs @@ -18,7 +18,7 @@ fn test_radau5_robertson_small_h() { params.set_tolerances(1e-2, 1e-2, None).unwrap(); // allocate the solver - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); // solve the ODE system let res = solver.solve(&mut y0, x0, x1, None, &mut args); diff --git a/russell_ode/tests/test_radau5_van_der_pol.rs b/russell_ode/tests/test_radau5_van_der_pol.rs index d9ef7634..aa09396d 100644 --- a/russell_ode/tests/test_radau5_van_der_pol.rs +++ b/russell_ode/tests/test_radau5_van_der_pol.rs @@ -12,7 +12,7 @@ fn test_radau5_van_der_pol() { params.step.h_ini = 1e-6; // allocate the solver - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); // enable dense output solver diff --git a/russell_ode/tests/test_radau5_van_der_pol_debug.rs b/russell_ode/tests/test_radau5_van_der_pol_debug.rs index 1e1d02c1..3b98fdbe 100644 --- a/russell_ode/tests/test_radau5_van_der_pol_debug.rs +++ b/russell_ode/tests/test_radau5_van_der_pol_debug.rs @@ -13,7 +13,7 @@ fn test_radau5_van_der_pol_debug() { params.debug = true; // allocate the solver - let mut solver = OdeSolver::new(params, &system).unwrap(); + let mut solver = OdeSolver::new(params, system).unwrap(); // solve the ODE system solver.solve(&mut y0, x0, x1, None, &mut args).unwrap();