[wip] Store Output in OdeSolver

.vscode/settings.json

@@ -151,5 +151,10 @@
     "zscal",
     "zsyrk",
     "zticks"
+  ],
+  "rust-analyzer.cargo.features": [
+    "intel_mkl",
+    "local_suitesparse",
+    "with_mumps"
   ]
 }

README.md

@@ -600,10 +600,9 @@ fn main() -> Result<(), StrError> {
     let mut solver = OdeSolver::new(params, &system)?;
 
     // enable dense output
-    let mut out = Output::new();
     let h_out = 0.01;
     let selected_y_components = &[0, 1];
-    out.set_dense_recording(true, h_out, selected_y_components)?;
+    solver.enable_output().set_dense_recording(true, h_out, selected_y_components)?;
 
     // solve the problem
     solver.solve(&mut y0, x0, x1, None, Some(&mut out), &mut args)?;

russell_ode/examples/amplifier1t_radau5.rs

@@ -24,14 +24,15 @@ fn main() -> Result<(), StrError> {
     let mut solver = OdeSolver::new(params, &system)?;
 
     // enable dense output
-    let mut out = Output::new();
     let h_out = 0.0001;
     let selected_y_components = &[0, 4];
-    out.set_dense_recording(true, h_out, selected_y_components)?;
+    solver
+        .enable_output()
+        .set_dense_recording(true, h_out, selected_y_components)?;
 
     // solve the problem
     let x1 = 0.2;
-    solver.solve(&mut y0, x0, x1, None, Some(&mut out), &mut args)?;
+    solver.solve(&mut y0, x0, x1, None, &mut args)?;
 
     // print the results and stats
     let y_ref = &[
@@ -74,9 +75,9 @@ fn main() -> Result<(), StrError> {
         .set_marker_style("+")
         .set_line_style("None");
 
-    curve1.draw(&out.dense_x, out.dense_y.get(&0).unwrap());
+    curve1.draw(&solver.out().dense_x, solver.out().dense_y.get(&0).unwrap());
     curve2.draw(&math.x, &math.y0);
-    curve3.draw(&out.dense_x, out.dense_y.get(&4).unwrap());
+    curve3.draw(&solver.out().dense_x, solver.out().dense_y.get(&4).unwrap());
     curve4.draw(&math.x, &math.y4);
 
     // save figure

russell_ode/examples/arenstorf_dopri8.rs

@@ -19,19 +19,22 @@ fn main() -> Result<(), StrError> {
     // get the ODE system
     let (system, x0, mut y0, x1, mut args, y_ref) = Samples::arenstorf();
 
-    // solver
+    // set configuration parameters
     let params = Params::new(Method::DoPri8);
+
+    // allocate the solver
     let mut solver = OdeSolver::new(params, &system)?;
 
     // enable dense output
-    let mut out = Output::new();
     let h_out = 0.01;
     let selected_y_components = &[0, 1];
-    out.set_dense_recording(true, h_out, selected_y_components)?;
+    solver
+        .enable_output()
+        .set_dense_recording(true, h_out, selected_y_components)?;
 
     // solve the problem
     let y = &mut y0;
-    solver.solve(y, x0, x1, None, Some(&mut out), &mut args)?;
+    solver.solve(y, x0, x1, None, &mut args)?;
 
     // print the results and stats
     println!("y_russell     = {:?}", y.as_data());
@@ -46,7 +49,10 @@ fn main() -> Result<(), StrError> {
     let mut curve2 = Curve::new();
     curve1.set_label("russell");
     curve2.set_label("mathematica");
-    curve1.draw(out.dense_y.get(&0).unwrap(), out.dense_y.get(&1).unwrap());
+    curve1.draw(
+        solver.out().dense_y.get(&0).unwrap(),
+        solver.out().dense_y.get(&1).unwrap(),
+    );
     curve2.set_marker_style(".").set_line_style("None");
     curve2.draw(&math.y0, &math.y1);
 

russell_ode/examples/brusselator_ode_dopri8.rs

@@ -22,18 +22,21 @@ fn main() -> Result<(), StrError> {
     // final x
     let x1 = 20.0;
 
-    // solver
+    // set configuration parameters
     let params = Params::new(Method::DoPri8);
+
+    // allocate the solver
     let mut solver = OdeSolver::new(params, &system)?;
 
     // enable dense output
-    let mut out = Output::new();
     let h_out = 0.01;
     let selected_y_components = &[0, 1];
-    out.set_dense_recording(true, h_out, selected_y_components)?;
+    solver
+        .enable_output()
+        .set_dense_recording(true, h_out, selected_y_components)?;
 
     // solve the problem
-    solver.solve(&mut y0, x0, x1, None, Some(&mut out), &mut args)?;
+    solver.solve(&mut y0, x0, x1, None, &mut args)?;
 
     // print the results and stats
     println!("y_russell     = {:?}", y0.as_data());
@@ -48,7 +51,10 @@ fn main() -> Result<(), StrError> {
     let mut curve2 = Curve::new();
     curve1.set_label("russell");
     curve2.set_label("mathematica");
-    curve1.draw(out.dense_y.get(&0).unwrap(), out.dense_y.get(&1).unwrap());
+    curve1.draw(
+        solver.out().dense_y.get(&0).unwrap(),
+        solver.out().dense_y.get(&1).unwrap(),
+    );
     curve2.set_marker_style(".").set_line_style("None");
     curve2.draw(&math.y0, &math.y1);
 

russell_ode/examples/brusselator_ode_fix_step.rs

@@ -55,7 +55,7 @@ fn main() -> Result<(), StrError> {
         print!("{:>w$}", name, w = w1);
         for i in 0..hh.len() {
             let mut y = y0.clone();
-            solver.solve(&mut y, x0, x1, Some(hh[i]), None, &mut args).unwrap();
+            solver.solve(&mut y, x0, x1, Some(hh[i]), &mut args).unwrap();
 
             // compare with the reference solution
             let (_, err) = vec_max_abs_diff(&y, &y_ref)?;

russell_ode/examples/brusselator_ode_var_step.rs

@@ -61,7 +61,7 @@ fn main() -> Result<(), StrError> {
 
             // call solve
             let mut y = y0.clone();
-            solver.solve(&mut y, x0, x1, None, None, &mut args).unwrap();
+            solver.solve(&mut y, x0, x1, None, &mut args).unwrap();
 
             // compare with the reference solution
             let (_, err) = vec_max_abs_diff(&y, &y_ref)?;

russell_ode/examples/brusselator_pde_2nd_comparison.rs

@@ -21,10 +21,12 @@ fn main() {
     let mut params = Params::new(Method::Radau5);
     params.set_tolerances(1e-4, 1e-4, None).unwrap();
 
-    // solve the ODE system
+    // allocate the solver
     let mut solver = OdeSolver::new(params, &system).unwrap();
+
+    // solve the ODE system
     let mut yy = yy0.clone();
-    solver.solve(&mut yy, t0, t1, None, None, &mut args).unwrap();
+    solver.solve(&mut yy, t0, t1, None, &mut args).unwrap();
 
     // get statistics
     let stat = solver.stats();

russell_ode/examples/brusselator_pde_radau5.rs

@@ -25,14 +25,15 @@ fn main() -> Result<(), StrError> {
     let mut params = Params::new(Method::Radau5);
     params.set_tolerances(1e-4, 1e-4, None)?;
 
+    // allocate the solver
+    let mut solver = OdeSolver::new(params, &system)?;
+
     // output
-    let mut out = Output::new();
     let h_out = 0.5;
-    out.set_dense_file_writing(true, h_out, PATH_KEY)?;
+    solver.enable_output().set_dense_file_writing(true, h_out, PATH_KEY)?;
 
     // solve the ODE system
-    let mut solver = OdeSolver::new(params, &system)?;
-    solver.solve(&mut yy0, t0, t1, None, Some(&mut out), &mut args)?;
+    solver.solve(&mut yy0, t0, t1, None, &mut args)?;
 
     // get statistics
     let stat = solver.stats();

russell_ode/examples/brusselator_pde_radau5_2nd.rs

@@ -26,14 +26,15 @@ fn main() -> Result<(), StrError> {
     let mut params = Params::new(Method::Radau5);
     params.set_tolerances(1e-4, 1e-4, None)?;
 
+    // allocate the solver
+    let mut solver = OdeSolver::new(params, &system)?;
+
     // output
-    let mut out = Output::new();
     let h_out = 1.0;
-    out.set_dense_file_writing(true, h_out, PATH_KEY)?;
+    solver.enable_output().set_dense_file_writing(true, h_out, PATH_KEY)?;
 
     // solve the ODE system
-    let mut solver = OdeSolver::new(params, &system)?;
-    solver.solve(&mut yy0, t0, t1, None, Some(&mut out), &mut args)?;
+    solver.solve(&mut yy0, t0, t1, None, &mut args)?;
 
     // get statistics
     let stat = solver.stats();

russell_ode/examples/hairer_wanner_eq1.rs

@@ -30,25 +30,26 @@ fn main() -> Result<(), StrError> {
     let mut fweuler = OdeSolver::new(Params::new(Method::FwEuler), &system)?;
 
     // solve the problem with BwEuler and h = 0.5
-    let mut out1 = Output::new();
-    out1.set_step_recording(true, &[0]);
+    bweuler.enable_output().set_step_recording(true, &[0]);
     let h = 0.5;
     let mut y = y0.clone();
-    bweuler.solve(&mut y, x0, x1, Some(h), Some(&mut out1), &mut args)?;
+    bweuler.solve(&mut y, x0, x1, Some(h), &mut args)?;
 
     // solve the problem with FwEuler and h = 1.974/50.0
-    let mut out2 = Output::new();
-    out2.set_step_recording(true, &[0]);
+    fweuler.enable_output().set_step_recording(true, &[0]);
     let h = 1.974 / 50.0;
     let mut y = y0.clone();
-    fweuler.solve(&mut y, x0, x1, Some(h), Some(&mut out2), &mut args)?;
+    fweuler.solve(&mut y, x0, x1, Some(h), &mut args)?;
+
+    // save the results for later
+    let out2_x = fweuler.out().step_x.clone();
+    let out2_y = fweuler.out().step_y.get(&0).unwrap().clone();
 
     // solve the problem with FwEuler and h = 1.875/50.0
-    let mut out3 = Output::new();
-    out3.set_step_recording(true, &[0]);
+    fweuler.enable_output().clear().set_step_recording(true, &[0]);
     let h = 1.875 / 50.0;
     let mut y = y0.clone();
-    fweuler.solve(&mut y, x0, x1, Some(h), Some(&mut out3), &mut args)?;
+    fweuler.solve(&mut y, x0, x1, Some(h), &mut args)?;
 
     // analytical solution
     let mut y_aux = Vector::new(system.get_ndim());
@@ -67,17 +68,15 @@ fn main() -> Result<(), StrError> {
     let mut curve1 = Curve::new();
     curve1
         .set_label("BwEuler h = 0.5")
-        .draw(&out1.step_x, out1.step_y.get(&0).unwrap());
+        .draw(&bweuler.out().step_x, bweuler.out().step_y.get(&0).unwrap());
 
     // FwEuler curves
     let mut curve2 = Curve::new();
     let mut curve3 = Curve::new();
-    curve2
-        .set_label("FwEuler h = 1.974/50")
-        .draw(&out2.step_x, out2.step_y.get(&0).unwrap());
+    curve2.set_label("FwEuler h = 1.974/50").draw(&out2_x, &out2_y);
     curve3
         .set_label("FwEuler h = 1.875/50")
-        .draw(&out3.step_x, out3.step_y.get(&0).unwrap());
+        .draw(&fweuler.out().step_x, fweuler.out().step_y.get(&0).unwrap());
 
     // save figure
     let mut plot = Plot::new();

russell_ode/examples/pde_1d_heat_spectral_collocation.rs

@@ -106,7 +106,7 @@ fn run(
         |x| f64::sin((x + 1.0) * PI));
 
     // solve the problem
-    ode.solve(&mut uu, t0, t1, None, None, &mut args)?;
+    ode.solve(&mut uu, t0, t1, None,  &mut args)?;
 
     // print stats
     if print_stats {

russell_ode/examples/robertson.rs

@@ -50,27 +50,28 @@ fn main() -> Result<(), StrError> {
     let sel = 1;
 
     // solve the problem with Radau5
-    let mut out1 = Output::new();
-    out1.set_step_recording(true, &[sel]);
+    radau5.enable_output().set_step_recording(true, &[sel]);
     let mut y = y0.clone();
-    radau5.solve(&mut y, x0, x1, None, Some(&mut out1), &mut args)?;
+    radau5.solve(&mut y, x0, x1, None, &mut args)?;
     println!("{}", radau5.stats());
     let n_accepted1 = radau5.stats().n_accepted;
 
     // solve the problem with DoPri5 and Tol = 1e-2
-    let mut out2 = Output::new();
-    out2.set_step_recording(true, &[sel]);
+    dopri5.enable_output().set_step_recording(true, &[sel]);
     let mut y = y0.clone();
-    dopri5.solve(&mut y, x0, x1, None, Some(&mut out2), &mut args)?;
+    dopri5.solve(&mut y, x0, x1, None, &mut args)?;
     println!("\nTol = 1e-2\n{}", dopri5.stats());
     let n_accepted2 = dopri5.stats().n_accepted;
 
-    // solve the problem with DoPri5 and Tol = 1e-3
-    let mut out3 = Output::new();
-    out3.set_step_recording(true, &[sel]);
+    // save the results for later
+    let out2_x = dopri5.out().step_x.clone();
+    let out2_y = dopri5.out().step_y.get(&sel).unwrap().clone();
+
+    // solve the problem again with DoPri5 and Tol = 1e-3
+    dopri5.enable_output().clear().set_step_recording(true, &[sel]);
     let mut y = y0.clone();
     dopri5.update_params(params3)?;
-    dopri5.solve(&mut y, x0, x1, None, Some(&mut out3), &mut args)?;
+    dopri5.solve(&mut y, x0, x1, None, &mut args)?;
     println!("\nTol = 1e-3\n{}", dopri5.stats());
     let n_accepted3 = dopri5.stats().n_accepted;
 
@@ -79,19 +80,19 @@ fn main() -> Result<(), StrError> {
     curve1
         .set_label(&format!("Radau5, n_accepted = {}", n_accepted1))
         .set_marker_style("o")
-        .draw(&out1.step_x, out1.step_y.get(&sel).unwrap());
+        .draw(&radau5.out().step_x, radau5.out().step_y.get(&sel).unwrap());
 
     // DoPri5 curves
     let mut curve2 = Curve::new();
     let mut curve3 = Curve::new();
     let mut curve4 = Curve::new();
     curve2
         .set_label(&format!("DoPri5, Tol = 1e-2, n_accepted = {}", n_accepted2))
-        .draw(&out2.step_x, out2.step_y.get(&sel).unwrap());
+        .draw(&out2_x, &out2_y);
     curve3
         .set_label(&format!("DoPri5, Tol = 1e-3, n_accepted = {}", n_accepted3))
-        .draw(&out3.step_x, out3.step_y.get(&sel).unwrap());
-    curve4.draw(&out3.step_x, &out3.step_h);
+        .draw(&dopri5.out().step_x, dopri5.out().step_y.get(&sel).unwrap());
+    curve4.draw(&dopri5.out().step_x, &dopri5.out().step_h);
 
     // save figures
     let mut plot1 = Plot::new();

russell_ode/examples/simple_ode_single_equation.rs

@@ -27,7 +27,7 @@ fn main() -> Result<(), StrError> {
     // solve from x = 0 to x = 1
     let x1 = 1.0;
     let mut args = 0;
-    solver.solve(&mut y, x, x1, None, None, &mut args)?;
+    solver.solve(&mut y, x, x1, None, &mut args)?;
     println!("y =\n{}", y);
 
     // check the results

russell_ode/examples/simple_system_with_mass.rs

@@ -47,7 +47,7 @@ fn main() -> Result<(), StrError> {
     // solve from x = 0 to x = 20
     let x1 = 20.0;
     let mut args = 0;
-    solver.solve(&mut y, x, x1, None, None, &mut args)?;
+    solver.solve(&mut y, x, x1, None, &mut args)?;
     println!("y =\n{}", y);
 
     // check the results