Ode make system copy

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;

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;
@@ -259,16 +259,17 @@ 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);
-    let mut solver = OdeSolver::new(params, &system)?;
+    let mut solver = OdeSolver::new(params, system)?;
 
     // initial values
     let x = 0.0;

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;

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;

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;

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()];

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()];

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();

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;

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;

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

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();

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;

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;

russell_ode/examples/simple_system_with_mass.rs

@@ -30,16 +30,17 @@ 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);
-    let mut solver = OdeSolver::new(params, &system)?;
+    let mut solver = OdeSolver::new(params, system)?;
 
     // initial values
     let x = 0.0;

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;

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];

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;

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);

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,22 +34,23 @@ 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
         } else {
             system.jac_nnz
         };
         let nnz = jac_nnz + ndim; // +ndim corresponds to the diagonal I matrix
+        let sym = system.symmetric;
         EulerBackward {
             params,
             system,
             k: Vector::new(ndim),
             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<A> 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]);

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]);