Merge pull request #2262 from oscardssmith/os/fix-inf-norm-for-in-place

fix `eigen_est` for in place explicit solvers with `Complex` `u0`
SciML · Jun 27, 2024 · e17f08f · e17f08f
2 parents 2f1f172 + 62c5a7f
commit e17f08f
Show file tree

Hide file tree

Showing 4 changed files with 11 additions and 7 deletions.
diff --git a/src/perform_step/explicit_rk_perform_step.jl b/src/perform_step/explicit_rk_perform_step.jl
@@ -222,7 +222,7 @@ end
     if integrator.alg isa CompositeAlgorithm
         # Hairer II, page 22 modified to use Inf norm
         @.. broadcast=false utilde=abs((kk[end] - kk[end - 1]) / (u - tmp))
-        integrator.eigen_est = integrator.opts.internalnorm(maximum(utilde), t)
+        integrator.eigen_est = integrator.opts.internalnorm(norm(utilde, Inf), t)
     end
 
     if integrator.opts.adaptive

diff --git a/src/perform_step/low_order_rk_perform_step.jl b/src/perform_step/low_order_rk_perform_step.jl
@@ -838,7 +838,7 @@ end
         g6 = tmp
         # Hairer II, page 22 modified to use Inf norm
         @.. broadcast=false thread=thread utilde=abs((k7 - k6) / (g7 - g6))
-        integrator.eigen_est = integrator.opts.internalnorm(maximum(utilde), t)
+        integrator.eigen_est = integrator.opts.internalnorm(norm(utilde, Inf), t)
     end
     if integrator.opts.adaptive
         @.. broadcast=false thread=thread utilde=dt * (btilde1 * k1 + btilde2 * k2 +
@@ -956,7 +956,7 @@ end
         g7 = u
         # Hairer II, page 22 modified to use Inf norm
         @.. broadcast=false thread=thread utilde=abs((k7 - k6) / (g7 - g6))
-        integrator.eigen_est = integrator.opts.internalnorm(maximum(utilde) * oneunit(t), t)
+        integrator.eigen_est = integrator.opts.internalnorm(norm(utilde, Inf) * oneunit(t), t)
     end
     if integrator.opts.adaptive
         @.. broadcast=false thread=thread utilde=dt * (btilde1 * k1 + btilde3 * k3 +

diff --git a/src/perform_step/verner_rk_perform_step.jl b/src/perform_step/verner_rk_perform_step.jl
@@ -163,7 +163,7 @@ end
         g9 = u
         g8 = tmp
         @.. broadcast=false thread=thread rtmp=abs((k9 - k8) / (g9 - g8))
-        integrator.eigen_est = integrator.opts.internalnorm(maximum(rtmp), t)
+        integrator.eigen_est = integrator.opts.internalnorm(norm(rtmp, Inf), t)
     end
     if integrator.opts.adaptive
         @.. broadcast=false thread=thread utilde=dt * (btilde1 * k1 + btilde4 * k4 +
@@ -411,7 +411,7 @@ end
         g10 = u
         g9 = tmp
         @.. broadcast=false thread=thread rtmp=abs((k10 - k9) / (g10 - g9))
-        integrator.eigen_est = integrator.opts.internalnorm(maximum(rtmp), t)
+        integrator.eigen_est = integrator.opts.internalnorm(norm(rtmp, Inf), t)
     end
     if integrator.opts.adaptive
         @.. broadcast=false thread=thread utilde=dt * (btilde1 * k1 + btilde4 * k4 +
@@ -731,7 +731,7 @@ end
         g13 = u
         g12 = tmp
         @.. broadcast=false thread=thread rtmp=abs((k13 - k12) / (g13 - g12))
-        integrator.eigen_est = integrator.opts.internalnorm(maximum(rtmp), t)
+        integrator.eigen_est = integrator.opts.internalnorm(norm(rtmp, Inf), t)
     end
     @.. broadcast=false thread=thread u=uprev +
                                         dt *
@@ -1136,7 +1136,7 @@ end
         g16 = u
         g15 = tmp
         @.. broadcast=false thread=thread rtmp=abs((k16 - k15) / (g16 - g15))
-        integrator.eigen_est = integrator.opts.internalnorm(maximum(rtmp), t)
+        integrator.eigen_est = integrator.opts.internalnorm(norm(rtmp, Inf), t)
     end
     @.. broadcast=false thread=thread u=uprev +
                                         dt *

diff --git a/test/interface/composite_algorithm_test.jl b/test/interface/composite_algorithm_test.jl
@@ -72,3 +72,7 @@ sol = solve(prob, alg = AutoVern7(Rodas5()))
 sol = solve(prob,
     alg = OrdinaryDiffEq.AutoAlgSwitch(ExplicitRK(constructVerner7()), Rodas5()))
 @test sol.t[end] == 1000.0
+
+prob = remake(prob_ode_2Dlinear, u0=rand(ComplexF64, 2, 2))
+sol = solve(prob, AutoTsit5(Rosenbrock23(autodiff=false))) # Complex and AD don't mix
+@test sol.retcode == ReturnCode.Success