❗ Panic when including mass matrix in v0.6.0

[MasoudAE21]

Hi all,

I'm encountering a panic when running a test that uses the v0.6.0, specifically when I include a mass matrix in the problem setup.

The error is:

thread 'tests::test_diffsol' panicked at ...
called `Option::unwrap()` on a `None` value
--> diffsol-0.6.0/src/matrix/dense_nalgebra_serial.rs:178:59

This happens only when I define a mass matrix function — if I omit that and solve a standard ODE, everything works fine.

---

🧪 Minimal Test Case

Here's the test function that triggers the panic:

fn test_diffsol() {
    let rhs_funcion = |x: &V, p: &V, t: T, y: &mut V| {
        y[0] = 2.0;
        println!("RHS: x: {}, y: {}, t: {}", x[0], y[0], t);
    };

    let jacobian_function = |x: &V, p: &V, t: T, v: &V, y: &mut V| {
        y[0] = 0.0;
        println!("JAC: y: {}; x: {}; t: {}, v: {}", y[0], x[0], t, v[0]);
    };

    let init_func = |p: &V, t: T, y: &mut V| {
        y[0] = 0.0;
        println!("INIT: y: {}; t: {}", y[0], t);
    };

    let mass_func = |v: &V, p: &V, t: T, beta: T, y: &mut V| {
        y[0] = v[0];
        println!("MASS: x: {}, y: {}, t: {}, beta: {}", v[0], y[0], t, beta);
    };

    let options = SemiExplicitDAESolverOptions {
        rel_tol: 1e-10,
        abs_tol: 1e-10,
        max_steps: 10000,
        min_step: 1e-6,
    };

    let problem = OdeBuilder::<M>::new()
        .rhs_implicit(rhs_funcion, jacobian_function)
        .mass(mass_func)
        .init(init_func, 1)
        .build()
        .unwrap();

    let mut diffsol_solver = problem.esdirk34::<LS>().unwrap_or_else(|e| {
        panic!("Error creating esdirk34 solver: {}", e);
    });

    let (ans, t) = diffsol_solver.solve(10.0).unwrap();

    println!("{:#?}", t);
}

---

🧠 What I Know

• The panic occurs inside dense_nalgebra_serial.rs at line 178, which seems to relate to gather() logic in matrix handling.
• If I remove the .mass(mass_func) call, the solver works fine for basic ODEs.
• Stack trace shows the issue arises during Matrix::split() inside solver initialization.

---

🙏 Questions

1. Is there a specific requirement for matrix shape or sparsity when defining the mass function?
2. Could I be missing a call that informs the solver about structure or dimensions?
3. Should .mass() only be used with certain solvers or problem types?

Any advice or clarification would be much appreciated — including if there's a known workaround or internal check I can add before calling .solve().

Thanks in advance!

[martinjrobins]

Thanks for reporting this @MasoudAE21, its a bug due to having a mass matrix with no algebraic variables, I've put in a fix here: #153

1. Is there a specific requirement for matrix shape or sparsity when defining the mass function?

The mass matrix needs to be a n x n matrix where n is the number of states. Note that you don't provide the mass matrix directly, your function should assign into y the result of M * v + beta * y, which will be a vector of dim n

1. Could I be missing a call that informs the solver about structure or dimensions?

Nope, this was just a bug :)

1. Should .mass() only be used with certain solvers or problem types?

If your mass matrix is an identity then it is better not to use .mass()