fix non-identity mass matrix with default solver and test

Author: oscardssmith

src/alg_utils.jl

@@ -1081,3 +1081,5 @@ is_mass_matrix_alg(alg::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}) = false
 is_mass_matrix_alg(alg::CompositeAlgorithm) = all(is_mass_matrix_alg, alg.algs)
 is_mass_matrix_alg(alg::RosenbrockAlgorithm) = true
 is_mass_matrix_alg(alg::NewtonAlgorithm) = !isesdirk(alg)
+# hack for the default alg
+is_mass_matrix_alg(alg::CompositeAlgorithm{<:Any, <:Tuple{Tsit5, Vern7, Rosenbrock23, Rodas5P, FBDF, FBDF}}) = true

src/composite_algs.jl

@@ -173,7 +173,7 @@ function default_autoswitch(AS::AutoSwitchCache, integrator)
         AS.is_stiffalg, AS.current) ?
                AS.count < 0 ? 1 : AS.count + 1 :
                AS.count > 0 ? -1 : AS.count - 1
-    if (!AS.is_stiffalg && AS.count > AS.maxstiffstep)
+    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)

test/interface/default_solver_tests.jl

@@ -18,6 +18,11 @@ vernsol = solve(prob_ode_2Dlinear, Vern7(), reltol=1e-10)
 @test sol.stats.nf == vernsol.stats.nf
 @test all(isequal(2), sol.alg_choice)
 
+prob_ode_linear_fast = ODEProblem(ODEFunction(f_2dlinear, mass_matrix=2*I(2)), rand(2), (0.0, 1.0), 1.01)
+sol = solve(prob_ode_linear_fast)
+@test all(isequal(3), sol.alg_choice)
+# for some reason the timestepping here is different from regular Rosenbrock23 (including the initial timestep)
+
 function rober(u, p, t)
     y₁, y₂, y₃ = u
     k₁, k₂, k₃ = p
@@ -68,3 +73,12 @@ for n in (100, ) # 600 should be added but currently is broken for unknown reaso
     @test sol.alg_choice[1] == 1
     @test sol.alg_choice[end] == stiffalg
 end
+
+function swaplinear(u, p, t)
+    [u[2], u[1]].*p
+end
+swaplinearf = ODEFunction(swaplinear, mass_matrix=ones(2,2)-I(2))
+prob_swaplinear = ODEProblem(swaplinearf, rand(2), (0., 1.) 1.01)
+sol = solve(prob_swaplinear, reltol=1e-7) # reltol must be set to avoid running into a bug with Rosenbrock23
+@test all(isequal(4), sol.alg_choice)
+# for some reason the timestepping here is different from regular Rodas5P (including the initial timestep)