Skip to content

Commit

Permalink
Merge pull request #113 from cpmech/ode-improve-jacobian-setting
Browse files Browse the repository at this point in the history
Ode improve jacobian setting
  • Loading branch information
cpmech authored May 15, 2024
2 parents 0427220 + 00f9688 commit 57fba72
Show file tree
Hide file tree
Showing 15 changed files with 463 additions and 664 deletions.
38 changes: 15 additions & 23 deletions russell_ode/README.md
Original file line number Diff line number Diff line change
Expand Up @@ -127,17 +127,10 @@ use russell_ode::prelude::*;
fn main() -> Result<(), StrError> {
// system
let ndim = 1;
let system = System::new(
ndim,
|f, x, y, _args: &mut NoArgs| {
f[0] = x + y[0];
Ok(())
},
no_jacobian,
HasJacobian::No,
None,
None,
);
let system = System::new(ndim, |f, x, y, _args: &mut NoArgs| {
f[0] = x + y[0];
Ok(())
});

// solver
let params = Params::new(Method::DoPri8);
Expand Down Expand Up @@ -235,20 +228,21 @@ See the code [simple_system_with_mass.rs](https://github.com/cpmech/russell/tree
```rust
use russell_lab::{vec_approx_eq, StrError, Vector};
use russell_ode::prelude::*;
use russell_sparse::CooMatrix;
use russell_sparse::{CooMatrix, Sym};

fn main() -> Result<(), StrError> {
// ODE system
let ndim = 3;
let jac_nnz = 4;
let mut system = System::new(
ndim,
|f: &mut Vector, x: f64, y: &Vector, _args: &mut NoArgs| {
f[0] = -y[0] + y[1];
f[1] = y[0] + y[1];
f[2] = 1.0 / (1.0 + x);
Ok(())
},
let mut system = System::new(ndim, |f: &mut Vector, x: f64, y: &Vector, _args: &mut NoArgs| {
f[0] = -y[0] + y[1];
f[1] = y[0] + y[1];
f[2] = 1.0 / (1.0 + x);
Ok(())
});
system.set_jacobian(
Some(jac_nnz),
Sym::No,
move |jj: &mut CooMatrix, alpha: f64, _x: f64, _y: &Vector, _args: &mut NoArgs| {
jj.reset();
jj.put(0, 0, alpha * (-1.0)).unwrap();
Expand All @@ -257,9 +251,6 @@ fn main() -> Result<(), StrError> {
jj.put(1, 1, alpha * (1.0)).unwrap();
Ok(())
},
HasJacobian::Yes,
Some(jac_nnz),
None,
);

// mass matrix
Expand Down Expand Up @@ -293,6 +284,7 @@ fn main() -> Result<(), StrError> {
println!("{}", solver.stats());
Ok(())
}

```

The output looks like:
Expand Down
6 changes: 1 addition & 5 deletions russell_ode/examples/pde_1d_heat_spectral_collocation.rs
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@ use plotpy::{Curve, Plot, SuperTitleParams};
use russell_lab::algo::{InterpGrid, InterpLagrange, InterpParams};
use russell_lab::math::{GOLDEN_RATIO, PI};
use russell_lab::{mat_vec_mul, StrError, Vector};
use russell_ode::{no_jacobian, HasJacobian, Method, NoArgs, OdeSolver, Params, System};
use russell_ode::{Method, NoArgs, OdeSolver, Params, System};

const PATH_KEY: &str = "/tmp/russell_ode/pde_1d_heat_spectral_collocation";

Expand Down Expand Up @@ -80,10 +80,6 @@ fn run(
dudt[nn] = 0.0; // homogeneous BCs
Ok(())
},
no_jacobian,
HasJacobian::No,
None,
None,
);

// ODE solver
Expand Down
15 changes: 4 additions & 11 deletions russell_ode/examples/simple_ode_single_equation.rs
Original file line number Diff line number Diff line change
Expand Up @@ -4,17 +4,10 @@ use russell_ode::prelude::*;
fn main() -> Result<(), StrError> {
// ODE system
let ndim = 1;
let system = System::new(
ndim,
|f, x, y, _args: &mut NoArgs| {
f[0] = x + y[0];
Ok(())
},
no_jacobian,
HasJacobian::No,
None,
None,
);
let system = System::new(ndim, |f, x, y, _args: &mut NoArgs| {
f[0] = x + y[0];
Ok(())
});

// solver
let params = Params::new(Method::DoPri8);
Expand Down
22 changes: 10 additions & 12 deletions russell_ode/examples/simple_system_with_mass.rs
Original file line number Diff line number Diff line change
@@ -1,19 +1,20 @@
use russell_lab::{vec_approx_eq, StrError, Vector};
use russell_ode::prelude::*;
use russell_sparse::CooMatrix;
use russell_sparse::{CooMatrix, Sym};

fn main() -> Result<(), StrError> {
// ODE system
let ndim = 3;
let jac_nnz = 4;
let mut system = System::new(
ndim,
|f: &mut Vector, x: f64, y: &Vector, _args: &mut NoArgs| {
f[0] = -y[0] + y[1];
f[1] = y[0] + y[1];
f[2] = 1.0 / (1.0 + x);
Ok(())
},
let mut system = System::new(ndim, |f: &mut Vector, x: f64, y: &Vector, _args: &mut NoArgs| {
f[0] = -y[0] + y[1];
f[1] = y[0] + y[1];
f[2] = 1.0 / (1.0 + x);
Ok(())
});
system.set_jacobian(
Some(jac_nnz),
Sym::No,
move |jj: &mut CooMatrix, alpha: f64, _x: f64, _y: &Vector, _args: &mut NoArgs| {
jj.reset();
jj.put(0, 0, alpha * (-1.0)).unwrap();
Expand All @@ -22,9 +23,6 @@ fn main() -> Result<(), StrError> {
jj.put(1, 1, alpha * (1.0)).unwrap();
Ok(())
},
HasJacobian::Yes,
Some(jac_nnz),
None,
);

// mass matrix
Expand Down
11 changes: 0 additions & 11 deletions russell_ode/src/enums.rs
Original file line number Diff line number Diff line change
@@ -1,14 +1,3 @@
/// Indicates whether the analytical Jacobian is available or not
pub enum HasJacobian {
/// Indicates that the analytical Jacobian is available
Yes,

/// Indicates that the analytical Jacobian is not available
///
/// **Note:** The [crate::no_jacobian] function can be conveniently used with the `No` option.
No,
}

/// Holds information about the numerical method to solve (approximate) ODEs
#[derive(Clone, Copy, Debug)]
pub struct Information {
Expand Down
31 changes: 11 additions & 20 deletions russell_ode/src/erk_dense_out.rs
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,6 @@ use crate::StrError;
use crate::{Method, System};
use crate::{DORMAND_PRINCE_5_D, DORMAND_PRINCE_8_AD, DORMAND_PRINCE_8_CD, DORMAND_PRINCE_8_D};
use russell_lab::Vector;
use russell_sparse::CooMatrix;

/// Handles the dense output of explicit Runge-Kutta methods
pub(crate) struct ErkDenseOut {
Expand Down Expand Up @@ -58,9 +57,9 @@ impl ErkDenseOut {
}

/// Updates the data and returns the number of function evaluations
pub(crate) fn update<'a, F, J, A>(
pub(crate) fn update<'a, F, A>(
&mut self,
system: &System<F, J, A>,
system: &System<F, A>,
x: f64,
y: &Vector,
h: f64,
Expand All @@ -70,7 +69,6 @@ impl ErkDenseOut {
) -> Result<usize, StrError>
where
F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
J: Fn(&mut CooMatrix, f64, f64, &Vector, &mut A) -> Result<(), StrError>,
{
let mut n_function_eval = 0;

Expand Down Expand Up @@ -240,7 +238,7 @@ impl ErkDenseOut {
#[cfg(test)]
mod tests {
use super::ErkDenseOut;
use crate::{no_jacobian, HasJacobian, Method, System};
use crate::{Method, System};
use russell_lab::Vector;

#[test]
Expand Down Expand Up @@ -309,21 +307,14 @@ mod tests {
count: usize,
fail: usize,
}
let mut system = System::new(
1,
|_f: &mut Vector, _x: f64, _y: &Vector, args: &mut Args| {
args.count += 1;
if args.count == args.fail {
Err("STOP")
} else {
Ok(())
}
},
no_jacobian,
HasJacobian::No,
None,
None,
);
let mut system = System::new(1, |_f: &mut Vector, _x: f64, _y: &Vector, args: &mut Args| {
args.count += 1;
if args.count == args.fail {
Err("STOP")
} else {
Ok(())
}
});
let mut out = ErkDenseOut::new(Method::DoPri8, system.ndim).unwrap();
let mut args = Args { count: 0, fail: 1 };
let h = 0.1;
Expand Down
60 changes: 26 additions & 34 deletions russell_ode/src/euler_backward.rs
Original file line number Diff line number Diff line change
@@ -1,19 +1,18 @@
use crate::StrError;
use crate::{OdeSolverTrait, Params, System, Workspace};
use russell_lab::{vec_copy, vec_rms_scaled, vec_update, Vector};
use russell_sparse::{numerical_jacobian, CooMatrix, LinSolver, SparseMatrix};
use russell_sparse::{numerical_jacobian, LinSolver, SparseMatrix};

/// Implements the backward Euler (implicit) solver (implicit, order 1, unconditionally stable)
pub(crate) struct EulerBackward<'a, F, J, A>
pub(crate) struct EulerBackward<'a, F, A>
where
F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
J: Fn(&mut CooMatrix, f64, f64, &Vector, &mut A) -> Result<(), StrError>,
{
/// Holds the parameters
params: Params,

/// ODE system
system: &'a System<F, J, A>,
system: &'a System<'a, F, A>,

/// Vector holding the function evaluation
///
Expand All @@ -36,13 +35,12 @@ where
solver: LinSolver<'a>,
}

impl<'a, F, J, A> EulerBackward<'a, F, J, A>
impl<'a, F, A> EulerBackward<'a, F, A>
where
F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
J: Fn(&mut CooMatrix, f64, f64, &Vector, &mut A) -> Result<(), StrError>,
{
/// Allocates a new instance
pub fn new(params: Params, system: &'a System<F, J, A>) -> Self {
pub fn new(params: Params, system: &'a System<'a, F, A>) -> Self {
let ndim = system.ndim;
let jac_nnz = if params.newton.use_numerical_jacobian {
ndim * ndim
Expand All @@ -63,10 +61,9 @@ where
}
}

impl<'a, F, J, A> OdeSolverTrait<A> for EulerBackward<'a, F, J, A>
impl<'a, F, A> OdeSolverTrait<A> for EulerBackward<'a, F, A>
where
F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
J: Fn(&mut CooMatrix, f64, f64, &Vector, &mut A) -> Result<(), StrError>,
{
/// Enables dense output
fn enable_dense_output(&mut self) -> Result<(), StrError> {
Expand Down Expand Up @@ -115,13 +112,13 @@ where

// calculate J_new := h J
let kk = self.kk.get_coo_mut().unwrap();
if self.params.newton.use_numerical_jacobian || !self.system.jac_available {
if self.params.newton.use_numerical_jacobian || self.system.jacobian.is_none() {
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)?;
} else {
(self.system.jacobian)(kk, h, x_new, y_new, args)?;
(self.system.jacobian.as_ref().unwrap())(kk, h, x_new, y_new, args)?;
}

// add diagonal entries => calculate K = h J_new - I
Expand Down Expand Up @@ -190,8 +187,9 @@ where
#[cfg(test)]
mod tests {
use super::EulerBackward;
use crate::{HasJacobian, Method, OdeSolverTrait, Params, Samples, System, Workspace};
use crate::{Method, OdeSolverTrait, Params, Samples, System, Workspace};
use russell_lab::{array_approx_eq, Vector};
use russell_sparse::Sym;

// Mathematica code:
//
Expand Down Expand Up @@ -383,28 +381,22 @@ mod tests {
struct Args {
count_f: usize,
}
let system = System::new(
1,
|f, _, _, args: &mut Args| {
f[0] = 1.0;
args.count_f += 1;
if args.count_f == 1 {
Err("f: stop (count = 1)")
} else if args.count_f == 4 {
Err("f: stop (count = 4; num-jacobian)")
} else {
Ok(())
}
},
|jj, alpha, _x, _y, _args: &mut Args| {
jj.reset();
jj.put(0, 0, alpha * (0.0)).unwrap();
Err("jj: stop")
},
HasJacobian::Yes,
None,
None,
);
let mut system = System::new(1, |f, _, _, args: &mut Args| {
f[0] = 1.0;
args.count_f += 1;
if args.count_f == 1 {
Err("f: stop (count = 1)")
} else if args.count_f == 4 {
Err("f: stop (count = 4; num-jacobian)")
} else {
Ok(())
}
});
system.set_jacobian(None, Sym::No, |jj, alpha, _x, _y, _args: &mut Args| {
jj.reset();
jj.put(0, 0, alpha * (0.0)).unwrap();
Err("jj: stop")
});
let params = Params::new(Method::BwEuler);
let mut solver = EulerBackward::new(params, &system);
let mut work = Workspace::new(Method::BwEuler);
Expand Down
Loading

0 comments on commit 57fba72

Please sign in to comment.