From 1d8ad8502f31225be2658e0bea67c5ca4ebbb802 Mon Sep 17 00:00:00 2001
From: oscarddssmith <oscar.smith@juliacomputing.com>
Date: Tue, 28 May 2024 16:02:55 -0400
Subject: [PATCH] fix default solver for non-diagonal mass matrices

---
 src/composite_algs.jl                  | 8 ++++----
 test/interface/default_solver_tests.jl | 2 +-
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/composite_algs.jl b/src/composite_algs.jl
index 462fa083f3..ba954dc818 100644
--- a/src/composite_algs.jl
+++ b/src/composite_algs.jl
@@ -144,13 +144,13 @@ function nonstiffchoice(reltol)
     Int(x)
 end
 
-function stiffchoice(reltol, len)
+function stiffchoice(reltol, len, mass_matrix)
     x = if len > MEDIUMSIZE
         DefaultSolverChoice.KrylovFBDF
     elseif len > SMALLSIZE
         DefaultSolverChoice.FBDF
     else
-        if reltol < LOW_TOL
+        if reltol < LOW_TOL || !isdiag(mass_matrix)
             DefaultSolverChoice.Rodas5P
         else
             DefaultSolverChoice.Rosenbrock23
@@ -166,7 +166,7 @@ function default_autoswitch(AS::AutoSwitchCache, integrator)
     # Chooose the starting method
     if AS.current == 0
         choice = if AS.stiffalgfirst || integrator.f.mass_matrix != I
-            stiffchoice(reltol, len)
+            stiffchoice(reltol, len, integrator.f.mass_matrix)
         else
             nonstiffchoice(reltol)
         end
@@ -184,7 +184,7 @@ function default_autoswitch(AS::AutoSwitchCache, integrator)
     if integrator.f.mass_matrix != I || (!AS.is_stiffalg && AS.count > AS.maxstiffstep)
         integrator.dt = dt * AS.dtfac
         AS.is_stiffalg = true
-        AS.current = stiffchoice(reltol, len)
+        AS.current = stiffchoice(reltol, len, integrator.f.mass_matrix)
     elseif (AS.is_stiffalg && AS.count < -AS.maxnonstiffstep)
         integrator.dt = dt / AS.dtfac
         AS.is_stiffalg = false
diff --git a/test/interface/default_solver_tests.jl b/test/interface/default_solver_tests.jl
index 07aff4a3d2..39aa74e027 100644
--- a/test/interface/default_solver_tests.jl
+++ b/test/interface/default_solver_tests.jl
@@ -80,6 +80,6 @@ function swaplinear(u, p, t)
 end
 swaplinearf = ODEFunction(swaplinear, mass_matrix = ones(2, 2) - I(2))
 prob_swaplinear = ODEProblem(swaplinearf, rand(2), (0.0, 1.0), 1.01)
-sol = solve(prob_swaplinear, reltol = 1e-7) # reltol must be set to avoid running into a bug with Rosenbrock23
+sol = solve(prob_swaplinear)
 @test all(isequal(4), sol.alg_choice)
 # for some reason the timestepping here is different from regular Rodas5P (including the initial timestep)