Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Rewrite symmetry enum #79

Merged
merged 3 commits into from
Mar 9, 2024
Merged

Rewrite symmetry enum #79

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
b9a8254
[sparse] Rename Symmetry to Sym
cpmech Mar 9, 2024
36b3ce1
[sparse] Rename sym flag
cpmech Mar 9, 2024
88d63db
[IMPORTANT] Simplify the Sym enum
cpmech Mar 9, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion README.md
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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