Skip to content

Commit

Permalink
Merge pull request #79 from cpmech/rewrite-symmetry-enum
Browse files Browse the repository at this point in the history 
Rewrite symmetry enum
  • Loading branch information
@cpmech
cpmech authored Mar 9, 2024
2 parents 051f0f7 + 88d63db commit 8aa2cca
Show file tree
Hide file tree
Showing 46 changed files with 714 additions and 973 deletions.
2 changes: 1 addition & 1 deletion README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -294,7 +294,7 @@ fn main() -> Result<(), StrError> {
// . -1 -3 2 .
// . . 1 . .
// . 4 2 . 1
let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, None)?;
let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, Sym::No)?;
coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate
coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate
coo.put(1, 0, 3.0)?;
Expand Down
3 changes: 1 addition & 2 deletions russell_ode/src/euler_backward.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,15 +50,14 @@ where
system.jac_nnz
};
let nnz = jac_nnz + ndim; // +ndim corresponds to the diagonal I matrix
let symmetry = Some(system.jac_symmetry);
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, symmetry).unwrap(),
kk: SparseMatrix::new_coo(ndim, ndim, nnz, system.jac_sym).unwrap(),
solver: LinSolver::new(params.newton.genie).unwrap(),
}
}
Expand Down
162 changes: 57 additions & 105 deletions russell_ode/src/radau5.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -110,7 +110,6 @@ where
/// Allocates a new instance
pub fn new(params: Params, system: &'a System<F, J, A>) -> Self {
let ndim = system.ndim;
let symmetry = Some(system.jac_symmetry);
let mass_nnz = match system.mass_matrix.as_ref() {
Some(mass) => mass.get_info().2,
None => ndim,
Expand All @@ -125,9 +124,9 @@ where
Radau5 {
params,
system,
jj: SparseMatrix::new_coo(ndim, ndim, jac_nnz, symmetry).unwrap(),
kk_real: SparseMatrix::new_coo(ndim, ndim, nnz, symmetry).unwrap(),
kk_comp: ComplexSparseMatrix::new_coo(ndim, ndim, nnz, symmetry).unwrap(),
jj: SparseMatrix::new_coo(ndim, ndim, jac_nnz, system.jac_sym).unwrap(),
kk_real: SparseMatrix::new_coo(ndim, ndim, nnz, system.jac_sym).unwrap(),
kk_comp: ComplexSparseMatrix::new_coo(ndim, ndim, nnz, system.jac_sym).unwrap(),
solver_real: LinSolver::new(params.newton.genie).unwrap(),
solver_comp: ComplexLinSolver::new(params.newton.genie).unwrap(),
reuse_jacobian: false,
Expand Down Expand Up @@ -847,109 +846,62 @@ mod tests {
assert_eq!(work.bench.n_function, n_fcn_correct);
}

#[test]
fn radau5_works_mass_matrix() {
// problem
let (system, data, mut args) = Samples::simple_system_with_mass_matrix(false);
let yfx = data.y_analytical.unwrap();
let ndim = system.ndim;

// allocate structs
let params = Params::new(Method::Radau5);
let mut solver = Radau5::new(params, &system);
let mut work = Workspace::new(Method::Radau5);

// message
println!("{:>4}{:>10}{:>10}{:>10}", "step", "err_y0", "err_y1", "err_y2");

// numerical approximation
let h = 0.1;
let mut x = data.x0;
let mut y = data.y0.clone();
let mut y_ana = Vector::new(ndim);
for n in 0..5 {
// call step
solver.step(&mut work, x, &y, h, &mut args).unwrap();

// important: update n_accepted (must precede `accept`)
work.bench.n_accepted += 1;

// call accept
solver.accept(&mut work, &mut x, &mut y, h, &mut args).unwrap();

// important: save previous stepsize and relative error (must succeed `accept`)
work.h_prev = h;
work.rel_error_prev = f64::max(params.step.rel_error_prev_min, work.rel_error);

// check the results
yfx(&mut y_ana, x);
let err_y0 = f64::abs(y[0] - y_ana[0]);
let err_y1 = f64::abs(y[1] - y_ana[1]);
let err_y2 = f64::abs(y[2] - y_ana[2]);
println!(
"{:>4}{}{}{}",
n,
format_scientific(err_y0, 10, 2),
format_scientific(err_y1, 10, 2),
format_scientific(err_y2, 10, 2)
);
assert!(err_y0 < 1e-9);
assert!(err_y1 < 1e-9);
assert!(err_y2 < 1e-8);
}
}

#[test]
#[serial]
fn radau5_works_mass_matrix_symmetric_mumps() {
// problem
let (system, data, mut args) = Samples::simple_system_with_mass_matrix(true);
let yfx = data.y_analytical.unwrap();
let ndim = system.ndim;

// allocate structs
let mut params = Params::new(Method::Radau5);
params.newton.genie = Genie::Mumps;
let mut solver = Radau5::new(params, &system);
let mut work = Workspace::new(Method::Radau5);

// message
println!("{:>4}{:>10}{:>10}{:>10}", "step", "err_y0", "err_y1", "err_y2");

// numerical approximation
let h = 0.1;
let mut x = data.x0;
let mut y = data.y0.clone();
let mut y_ana = Vector::new(ndim);
for n in 0..5 {
// call step
solver.step(&mut work, x, &y, h, &mut args).unwrap();

// important: update n_accepted (must precede `accept`)
work.bench.n_accepted += 1;

// call accept
solver.accept(&mut work, &mut x, &mut y, h, &mut args).unwrap();

// important: save previous stepsize and relative error (must succeed `accept`)
work.h_prev = h;
work.rel_error_prev = f64::max(params.step.rel_error_prev_min, work.rel_error);

// check the results
yfx(&mut y_ana, x);
let err_y0 = f64::abs(y[0] - y_ana[0]);
let err_y1 = f64::abs(y[1] - y_ana[1]);
let err_y2 = f64::abs(y[2] - y_ana[2]);
println!(
"{:>4}{}{}{}",
n,
format_scientific(err_y0, 10, 2),
format_scientific(err_y1, 10, 2),
format_scientific(err_y2, 10, 2)
);
assert!(err_y0 < 1e-9);
assert!(err_y1 < 1e-9);
assert!(err_y2 < 1e-8);
fn radau5_works_mass_matrix() {
for symmetric in [true, false] {
for genie in [Genie::Umfpack, Genie::Mumps] {
// problem
let (system, data, mut args) = Samples::simple_system_with_mass_matrix(symmetric, genie);
let yfx = data.y_analytical.unwrap();
let ndim = system.ndim;

// allocate structs
let mut params = Params::new(Method::Radau5);
params.newton.genie = genie;
let mut solver = Radau5::new(params, &system);
let mut work = Workspace::new(Method::Radau5);

// message
println!("\nsymmetric = {:?} --- {:?}", symmetric, genie);
println!("{:>4}{:>10}{:>10}{:>10}", "step", "err_y0", "err_y1", "err_y2");

// numerical approximation
let h = 0.1;
let mut x = data.x0;
let mut y = data.y0.clone();
let mut y_ana = Vector::new(ndim);
for n in 0..4 {
// call step
solver.step(&mut work, x, &y, h, &mut args).unwrap();

// important: update n_accepted (must precede `accept`)
work.bench.n_accepted += 1;

// call accept
solver.accept(&mut work, &mut x, &mut y, h, &mut args).unwrap();

// important: save previous stepsize and relative error (must succeed `accept`)
work.h_prev = h;
work.rel_error_prev = f64::max(params.step.rel_error_prev_min, work.rel_error);

// check the results
yfx(&mut y_ana, x);
let err_y0 = f64::abs(y[0] - y_ana[0]);
let err_y1 = f64::abs(y[1] - y_ana[1]);
let err_y2 = f64::abs(y[2] - y_ana[2]);
println!(
"{:>4}{}{}{}",
n,
format_scientific(err_y0, 10, 2),
format_scientific(err_y1, 10, 2),
format_scientific(err_y2, 10, 2)
);
assert!(err_y0 < 1e-9);
assert!(err_y1 < 1e-9);
assert!(err_y2 < 1e-8);
}
}
}
}
}
59 changes: 25 additions & 34 deletions russell_ode/src/samples.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@ use crate::StrError;
use crate::{HasJacobian, NoArgs, System};
use russell_lab::math::PI;
use russell_lab::Vector;
use russell_sparse::{CooMatrix, Symmetry};
use russell_sparse::{CooMatrix, Genie, Sym};

/// Holds data corresponding to a sample ODE problem
pub struct SampleData {
Expand Down Expand Up @@ -139,7 +139,9 @@ impl Samples {
///
/// # Input
///
/// * `lower_triangle` -- Considers the symmetry of the Jacobian and Mass matrices, and generates a lower triangular representation
/// * `symmetric` -- considers the symmetry of the Jacobian and Mass matrices
/// * `genie` -- if symmetric, this information is required to decide on the lower-triangle/full
/// representation which is consistent with the linear solver employed
///
/// # Output
///
Expand All @@ -152,7 +154,8 @@ impl Samples {
/// * `data: SampleData` -- holds the initial values
/// * `args: NoArgs` -- is a placeholder variable with the arguments to F and J
pub fn simple_system_with_mass_matrix(
lower_triangle: bool,
symmetric: bool,
genie: Genie,
) -> (
System<
impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>,
Expand All @@ -162,12 +165,9 @@ impl Samples {
SampleData,
NoArgs,
) {
// selected symmetry option (for both Mass and Jacobian matrices)
let symmetry = if lower_triangle {
Some(Symmetry::new_general_lower())
} else {
None
};
// selected symmetric option (for both Mass and Jacobian matrices)
let sym = genie.symmetry(symmetric);
let triangular = sym.triangular();

// initial values
let x0 = 0.0;
Expand All @@ -176,7 +176,7 @@ impl Samples {

// ODE system
let ndim = 3;
let jac_nnz = if lower_triangle { 3 } else { 4 };
let jac_nnz = if triangular { 3 } else { 4 };
let mut system = System::new(
ndim,
move |f: &mut Vector, x: f64, y: &Vector, _args: &mut NoArgs| {
Expand All @@ -188,7 +188,7 @@ impl Samples {
move |jj: &mut CooMatrix, _x: f64, _y: &Vector, m: f64, _args: &mut NoArgs| {
jj.reset();
jj.put(0, 0, m * (-1.0)).unwrap();
if !lower_triangle {
if !triangular {
jj.put(0, 1, m * (1.0)).unwrap();
}
jj.put(1, 0, m * (1.0)).unwrap();
Expand All @@ -197,14 +197,14 @@ impl Samples {
},
HasJacobian::Yes,
Some(jac_nnz),
symmetry,
if sym == Sym::No { None } else { Some(sym) },
);

// mass matrix
let mass_nnz = if lower_triangle { 4 } else { 5 };
let mass_nnz = if triangular { 4 } else { 5 };
system.init_mass_matrix(mass_nnz).unwrap();
system.mass_put(0, 0, 1.0).unwrap();
if !lower_triangle {
if !triangular {
system.mass_put(0, 1, 1.0).unwrap();
}
system.mass_put(1, 0, 1.0).unwrap();
Expand Down Expand Up @@ -789,7 +789,7 @@ mod tests {
use super::{NoArgs, Samples};
use crate::StrError;
use russell_lab::{deriv_central5, mat_approx_eq, vec_approx_eq, Matrix, Vector};
use russell_sparse::{CooMatrix, Symmetry};
use russell_sparse::{CooMatrix, Sym};

fn numerical_jacobian<F>(ndim: usize, x0: f64, y0: Vector, function: F, multiplier: f64) -> Matrix
where
Expand Down Expand Up @@ -843,8 +843,7 @@ mod tests {
}

// compute the analytical Jacobian matrix
let symmetry = Some(system.jac_symmetry);
let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, symmetry).unwrap();
let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap();
(system.jacobian)(&mut jj, data.x0, &data.y0, multiplier, &mut args).unwrap();

// compute the numerical Jacobian matrix
Expand All @@ -871,8 +870,7 @@ mod tests {
}

// compute the analytical Jacobian matrix
let symmetry = Some(system.jac_symmetry);
let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, symmetry).unwrap();
let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap();
(system.jacobian)(&mut jj, data.x0, &data.y0, multiplier, &mut args).unwrap();

// compute the numerical Jacobian matrix
Expand All @@ -899,8 +897,7 @@ mod tests {
}

// compute the analytical Jacobian matrix
let symmetry = Some(system.jac_symmetry);
let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, symmetry).unwrap();
let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap();
(system.jacobian)(&mut jj, data.x0, &data.y0, multiplier, &mut args).unwrap();

// compute the numerical Jacobian matrix
Expand All @@ -927,8 +924,7 @@ mod tests {
}

// compute the analytical Jacobian matrix
let symmetry = Some(system.jac_symmetry);
let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, symmetry).unwrap();
let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap();
(system.jacobian)(&mut jj, data.x0, &data.y0, multiplier, &mut args).unwrap();

// compute the numerical Jacobian matrix
Expand All @@ -947,8 +943,7 @@ mod tests {
let (system, data, mut args) = Samples::robertson();

// compute the analytical Jacobian matrix
let symmetry = Some(system.jac_symmetry);
let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, symmetry).unwrap();
let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap();
(system.jacobian)(&mut jj, data.x0, &data.y0, multiplier, &mut args).unwrap();

// compute the numerical Jacobian matrix
Expand All @@ -967,8 +962,7 @@ mod tests {
let (system, data, mut args) = Samples::van_der_pol(None, false);

// compute the analytical Jacobian matrix
let symmetry = Some(system.jac_symmetry);
let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, symmetry).unwrap();
let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap();
(system.jacobian)(&mut jj, data.x0, &data.y0, multiplier, &mut args).unwrap();

// compute the numerical Jacobian matrix
Expand All @@ -987,8 +981,7 @@ mod tests {
let (system, data, mut args) = Samples::van_der_pol(None, true);

// compute the analytical Jacobian matrix
let symmetry = Some(system.jac_symmetry);
let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, symmetry).unwrap();
let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap();
(system.jacobian)(&mut jj, data.x0, &data.y0, multiplier, &mut args).unwrap();

// compute the numerical Jacobian matrix
Expand All @@ -1007,8 +1000,7 @@ mod tests {
let (system, data, mut args) = Samples::arenstorf();

// compute the analytical Jacobian matrix
let symmetry = Some(system.jac_symmetry);
let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, symmetry).unwrap();
let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap();
(system.jacobian)(&mut jj, data.x0, &data.y0, multiplier, &mut args).unwrap();

// compute the numerical Jacobian matrix
Expand All @@ -1027,8 +1019,7 @@ mod tests {
let (system, data, mut args) = Samples::amplifier1t();

// compute the analytical Jacobian matrix
let symmetry = Some(system.jac_symmetry);
let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, symmetry).unwrap();
let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap();
(system.jacobian)(&mut jj, data.x0, &data.y0, multiplier, &mut args).unwrap();

// compute the numerical Jacobian matrix
Expand All @@ -1045,6 +1036,6 @@ mod tests {
println!("{}", mass.as_dense());
let ndim = system.ndim;
let nnz_mass = 5 + 4;
assert_eq!(mass.get_info(), (ndim, ndim, nnz_mass, Symmetry::No));
assert_eq!(mass.get_info(), (ndim, ndim, nnz_mass, Sym::No));
}
}
Loading

0 comments on commit 8aa2cca

Please sign in to comment.