Skip to content

Commit

Permalink
Merge pull request #117 from cpmech/ode-convert-system-fn-to-boxed
Browse files Browse the repository at this point in the history
Convert the generic F parameter to a boxed function trait
  • Loading branch information
cpmech authored May 22, 2024
2 parents 3f33c06 + 7f0121d commit caa2c4c
Show file tree
Hide file tree
Showing 10 changed files with 49 additions and 144 deletions.
1 change: 1 addition & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -12,3 +12,4 @@ good.mtx
bad.mtx
good.smat
bad.smat
rustc-ice-*.txt
9 changes: 3 additions & 6 deletions russell_ode/src/erk_dense_out.rs
Original file line number Diff line number Diff line change
Expand Up @@ -57,19 +57,16 @@ impl ErkDenseOut {
}

/// Updates the data and returns the number of function evaluations
pub(crate) fn update<'a, F, A>(
pub(crate) fn update<'a, A>(
&mut self,
system: &System<F, A>,
system: &System<A>,
x: f64,
y: &Vector,
h: f64,
w: &Vector,
k: &Vec<Vector>,
args: &mut A,
) -> Result<usize, StrError>
where
F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
{
) -> Result<usize, StrError> {
let mut n_function_eval = 0;

if self.method == Method::DoPri5 {
Expand Down
19 changes: 5 additions & 14 deletions russell_ode/src/euler_backward.rs
Original file line number Diff line number Diff line change
Expand Up @@ -4,15 +4,12 @@ use russell_lab::{vec_copy, vec_rms_scaled, vec_update, Vector};
use russell_sparse::{numerical_jacobian, LinSolver, SparseMatrix};

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

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

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

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

impl<'a, F, A> OdeSolverTrait<A> for EulerBackward<'a, F, A>
where
F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
{
impl<'a, A> OdeSolverTrait<A> for EulerBackward<'a, A> {
/// Enables dense output
fn enable_dense_output(&mut self) -> Result<(), StrError> {
Err("dense output is not available for the BwEuler method")
Expand Down
19 changes: 5 additions & 14 deletions russell_ode/src/euler_forward.rs
Original file line number Diff line number Diff line change
Expand Up @@ -6,12 +6,9 @@ use russell_lab::{vec_add, vec_copy, Vector};
///
/// **Warning:** This method is interesting for didactic purposes only
/// and should not be used in production codes.
pub(crate) struct EulerForward<'a, F, A>
where
F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
{
pub(crate) struct EulerForward<'a, A> {
/// ODE system
system: &'a System<'a, F, A>,
system: &'a System<'a, A>,

/// Vector holding the function evaluation
///
Expand All @@ -22,12 +19,9 @@ where
w: Vector,
}

impl<'a, F, A> EulerForward<'a, F, A>
where
F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
{
impl<'a, A> EulerForward<'a, A> {
/// Allocates a new instance
pub fn new(system: &'a System<'a, F, A>) -> Self {
pub fn new(system: &'a System<'a, A>) -> Self {
let ndim = system.ndim;
EulerForward {
system,
Expand All @@ -37,10 +31,7 @@ where
}
}

impl<'a, F, A> OdeSolverTrait<A> for EulerForward<'a, F, A>
where
F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
{
impl<'a, A> OdeSolverTrait<A> for EulerForward<'a, A> {
/// Enables dense output
fn enable_dense_output(&mut self) -> Result<(), StrError> {
Err("dense output is not available for the FwEuler method")
Expand Down
19 changes: 5 additions & 14 deletions russell_ode/src/explicit_runge_kutta.rs
Original file line number Diff line number Diff line change
Expand Up @@ -20,15 +20,12 @@ use russell_lab::{format_fortran, vec_copy, vec_update, Matrix, Vector};
/// 2. E. Hairer, G. Wanner (2002) Solving Ordinary Differential Equations II.
/// Stiff and Differential-Algebraic Problems. Second Revised Edition.
/// Corrected 2nd printing 2002. Springer Series in Computational Mathematics, 614p
pub(crate) struct ExplicitRungeKutta<'a, F, A>
where
F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
{
pub(crate) struct ExplicitRungeKutta<'a, A> {
/// Holds the parameters
params: Params,

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

/// Information such as implicit, embedded, etc.
info: Information,
Expand Down Expand Up @@ -78,12 +75,9 @@ where
dense_out: Option<ErkDenseOut>,
}

impl<'a, F, A> ExplicitRungeKutta<'a, F, A>
where
F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
{
impl<'a, A> ExplicitRungeKutta<'a, A> {
/// Allocates a new instance
pub fn new(params: Params, system: &'a System<'a, F, A>) -> Result<Self, StrError> {
pub fn new(params: Params, system: &'a System<'a, A>) -> Result<Self, StrError> {
// Runge-Kutta coefficients
#[rustfmt::skip]
let (aa, bb, cc) = match params.method {
Expand Down Expand Up @@ -153,10 +147,7 @@ where
}
}

impl<'a, F, A> OdeSolverTrait<A> for ExplicitRungeKutta<'a, F, A>
where
F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
{
impl<'a, A> OdeSolverTrait<A> for ExplicitRungeKutta<'a, A> {
/// Enables dense output
fn enable_dense_output(&mut self) -> Result<(), StrError> {
self.dense_out = Some(ErkDenseOut::new(self.params.method, self.system.ndim)?);
Expand Down
9 changes: 3 additions & 6 deletions russell_ode/src/ode_solver.rs
Original file line number Diff line number Diff line change
Expand Up @@ -133,13 +133,10 @@ impl<'a, A> OdeSolver<'a, A> {
///
/// # Generics
///
/// The generic arguments here are:
///
/// * `F` -- function to compute the `f` vector: `(f: &mut Vector, x: f64, y: &Vector, args: &mut A)`
/// * `A` -- generic argument to assist in the `F` and Jacobian functions. It may be simply [crate::NoArgs] indicating that no arguments are needed.
pub fn new<F>(params: Params, system: &'a System<'a, F, A>) -> Result<Self, StrError>
/// * `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<Self, StrError>
where
F: 'a + Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
A: 'a,
{
if system.mass_matrix.is_some() && params.method != Method::Radau5 {
Expand Down
6 changes: 2 additions & 4 deletions russell_ode/src/output.rs
Original file line number Diff line number Diff line change
Expand Up @@ -38,10 +38,8 @@ pub struct OutCount {
///
/// # Generics
///
/// The generic arguments are:
///
/// * `A` -- Is auxiliary argument for the `F`, `J`, `YxFunction`, and `OutCallback` functions.
/// It may be simply [crate::NoArgs] indicating that no arguments are effectively used.
/// * `A` -- generic argument to assist in the f(x,y) and Jacobian functions.
/// It may be simply [crate::NoArgs] indicating that no arguments are needed.
pub struct Output<'a, A> {
/// Indicates whether the solver called initialize or not
initialized: bool,
Expand Down
19 changes: 5 additions & 14 deletions russell_ode/src/radau5.rs
Original file line number Diff line number Diff line change
Expand Up @@ -24,15 +24,12 @@ use std::thread;
/// 2. E. Hairer, G. Wanner (2002) Solving Ordinary Differential Equations II.
/// Stiff and Differential-Algebraic Problems. Second Revised Edition.
/// Corrected 2nd printing 2002. Springer Series in Computational Mathematics, 614p
pub(crate) struct Radau5<'a, F, A>
where
F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
{
pub(crate) struct Radau5<'a, A> {
/// Holds the parameters
params: Params,

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

/// Holds the Jacobian matrix. J = df/dy
jj: SparseMatrix,
Expand Down Expand Up @@ -118,12 +115,9 @@ where
dw12: ComplexVector, // packed (dw1, dw2)
}

impl<'a, F, A> Radau5<'a, F, A>
where
F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
{
impl<'a, A> Radau5<'a, A> {
/// Allocates a new instance
pub fn new(params: Params, system: &'a System<'a, F, A>) -> Self {
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,
Expand Down Expand Up @@ -331,10 +325,7 @@ where
}
}

impl<'a, F, A> OdeSolverTrait<A> for Radau5<'a, F, A>
where
F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
{
impl<'a, A> OdeSolverTrait<A> for Radau5<'a, A> {
/// Enables dense output
fn enable_dense_output(&mut self) -> Result<(), StrError> {
Ok(())
Expand Down
64 changes: 11 additions & 53 deletions russell_ode/src/samples.rs
Original file line number Diff line number Diff line change
@@ -1,4 +1,3 @@
use crate::StrError;
use crate::{NoArgs, PdeDiscreteLaplacian2d, Side, System};
use russell_lab::math::PI;
use russell_lab::Vector;
Expand Down Expand Up @@ -42,7 +41,7 @@ impl Samples {
/// * `args` -- is a placeholder variable with the arguments to F and J
/// * `y_fn_x` -- is a function to compute the analytical solution
pub fn simple_equation_constant<'a>() -> (
System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>,
System<'a, NoArgs>,
f64,
Vector,
NoArgs,
Expand Down Expand Up @@ -153,7 +152,7 @@ impl Samples {
symmetric: bool,
genie: Genie,
) -> (
System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>,
System<'a, NoArgs>,
f64,
Vector,
NoArgs,
Expand Down Expand Up @@ -257,13 +256,7 @@ impl Samples {
/// * Hairer E, Nørsett, SP, Wanner G (2008) Solving Ordinary Differential Equations I.
/// Non-stiff Problems. Second Revised Edition. Corrected 3rd printing 2008. Springer Series
/// in Computational Mathematics, 528p
pub fn brusselator_ode<'a>() -> (
System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>,
f64,
Vector,
NoArgs,
Vector,
) {
pub fn brusselator_ode<'a>() -> (System<'a, NoArgs>, f64, Vector, NoArgs, Vector) {
// system
let ndim = 2;
let jac_nnz = 4;
Expand Down Expand Up @@ -502,16 +495,7 @@ impl Samples {
npoint: usize,
second_book: bool,
ignore_diffusion: bool,
) -> (
System<
'a,
impl Fn(&mut Vector, f64, &Vector, &mut PdeDiscreteLaplacian2d) -> Result<(), StrError>,
PdeDiscreteLaplacian2d,
>,
f64,
Vector,
PdeDiscreteLaplacian2d,
) {
) -> (System<'a, PdeDiscreteLaplacian2d>, f64, Vector, PdeDiscreteLaplacian2d) {
// constants
let (kx, ky) = (alpha, alpha);
let (xmin, xmax) = (0.0, 1.0);
Expand Down Expand Up @@ -663,14 +647,7 @@ impl Samples {
/// * Hairer E, Nørsett, SP, Wanner G (2008) Solving Ordinary Differential Equations I.
/// Non-stiff Problems. Second Revised Edition. Corrected 3rd printing 2008. Springer Series
/// in Computational Mathematics, 528p
pub fn arenstorf<'a>() -> (
System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>,
f64,
Vector,
f64,
NoArgs,
Vector,
) {
pub fn arenstorf<'a>() -> (System<'a, NoArgs>, f64, Vector, f64, NoArgs, Vector) {
// constants
const MU: f64 = 0.012277471;
const MD: f64 = 1.0 - MU;
Expand Down Expand Up @@ -788,7 +765,7 @@ impl Samples {
/// Stiff and Differential-Algebraic Problems. Second Revised Edition.
/// Corrected 2nd printing 2002. Springer Series in Computational Mathematics, 614p
pub fn hairer_wanner_eq1<'a>() -> (
System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>,
System<'a, NoArgs>,
f64,
Vector,
NoArgs,
Expand Down Expand Up @@ -861,12 +838,7 @@ impl Samples {
/// * Hairer E, Wanner G (2002) Solving Ordinary Differential Equations II.
/// Stiff and Differential-Algebraic Problems. Second Revised Edition.
/// Corrected 2nd printing 2002. Springer Series in Computational Mathematics, 614p
pub fn robertson<'a>() -> (
System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>,
f64,
Vector,
NoArgs,
) {
pub fn robertson<'a>() -> (System<'a, NoArgs>, f64, Vector, NoArgs) {
// system
let ndim = 3;
let mut system = System::new(ndim, |f: &mut Vector, _x: f64, y: &Vector, _args: &mut NoArgs| {
Expand Down Expand Up @@ -942,16 +914,7 @@ impl Samples {
/// * Hairer E, Wanner G (2002) Solving Ordinary Differential Equations II.
/// Stiff and Differential-Algebraic Problems. Second Revised Edition.
/// Corrected 2nd printing 2002. Springer Series in Computational Mathematics, 614p
pub fn van_der_pol<'a>(
epsilon: f64,
stationary: bool,
) -> (
System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>,
f64,
Vector,
f64,
NoArgs,
) {
pub fn van_der_pol<'a>(epsilon: f64, stationary: bool) -> (System<'a, NoArgs>, f64, Vector, f64, NoArgs) {
// constants
let x0 = 0.0;
let mut y0 = Vector::from(&[2.0, -0.6]);
Expand Down Expand Up @@ -1071,12 +1034,7 @@ impl Samples {
/// * Hairer E, Wanner G (2002) Solving Ordinary Differential Equations II.
/// Stiff and Differential-Algebraic Problems. Second Revised Edition.
/// Corrected 2nd printing 2002. Springer Series in Computational Mathematics, 614p
pub fn amplifier1t<'a>() -> (
System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>,
f64,
Vector,
NoArgs,
) {
pub fn amplifier1t<'a>() -> (System<'a, NoArgs>, f64, Vector, NoArgs) {
// constants
const ALPHA: f64 = 0.99;
const GAMMA: f64 = 1.0 - ALPHA;
Expand Down Expand Up @@ -1174,7 +1132,7 @@ impl Samples {
/// * Kreyszig, E (2011) Advanced engineering mathematics; in collaboration with Kreyszig H,
/// Edward JN 10th ed 2011, Hoboken, New Jersey, Wiley
pub fn kreyszig_eq6_page902<'a>() -> (
System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>,
System<'a, NoArgs>,
f64,
Vector,
NoArgs,
Expand Down Expand Up @@ -1250,7 +1208,7 @@ impl Samples {
/// * Kreyszig, E (2011) Advanced engineering mathematics; in collaboration with Kreyszig H,
/// Edward JN 10th ed 2011, Hoboken, New Jersey, Wiley
pub fn kreyszig_ex4_page920<'a>() -> (
System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>,
System<'a, NoArgs>,
f64,
Vector,
NoArgs,
Expand Down
Loading

0 comments on commit caa2c4c

Please sign in to comment.