fixes

Author: Shreyas-Ekanathan

lib/OrdinaryDiffEqFIRK/src/firk_addsteps.jl

@@ -413,8 +413,8 @@ if integrator.EEst <= oneunit(integrator.EEst)
     if alg.extrapolant != :constant
         integrator.k[3] = (z2 - z3)/c2m1
         tmp = @.. (z1 - z2)/c1mc2
-        integrator.k[4] = (tmp - cache.cont1)/c1m1
-        integrator.k[5] = cache.cont2-(tmp - z1 / c1) / c2
+        integrator.k[4] = (tmp - integrator.k[3])/c1m1
+        integrator.k[5] = integrator.k[4]-(tmp - z1 / c1) / c2
     end
 end
 

lib/OrdinaryDiffEqFIRK/src/firk_interpolants.jl

@@ -7,6 +7,7 @@ FIRK_WITH_INTERPOLATIONS = Union{RadauIIA3ConstantCache, RadauIIA3Cache, RadauII
     @unpack cont1, cont2 = cache
     @unpack c1 = cache.tab
     c1m1 = c1 - 1
+    Θdt = 1 - Θ
     @.. y₁ - Θdt * (k[3] - (Θdt + c1m1) * k[4])
 end
 

lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl

@@ -63,7 +63,7 @@ function initialize!(integrator, cache::RadauIIA3Cache)
 end
 
 function initialize!(integrator, cache::RadauIIA5ConstantCache)
-    integrator.kshortsize = 2 
+    integrator.kshortsize = 5 
     integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
     integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
     OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
@@ -72,8 +72,9 @@ function initialize!(integrator, cache::RadauIIA5ConstantCache)
     integrator.fsallast = zero(integrator.fsalfirst)
     integrator.k[1] = integrator.fsalfirst
     integrator.k[2] = integrator.fsallast
-    integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t)
-    OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
+    integrator.k[3] = zero(integrator.fsalfirst)
+    integrator.k[4] = zero(integrator.fsalfirst)
+    integrator.k[5] = zero(integrator.fsalfirst)
     nothing
 end
 
@@ -110,8 +111,11 @@ function initialize!(integrator, cache::RadauIIA9ConstantCache)
     integrator.fsallast = zero(integrator.fsalfirst)
     integrator.k[1] = integrator.fsalfirst
     integrator.k[2] = integrator.fsallast
-    integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t)
-    OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
+    integrator.k[3] = zero(integrator.fsalfirst)
+    integrator.k[4] = zero(integrator.fsalfirst)
+    integrator.k[5] = zero(integrator.fsalfirst)
+    integrator.k[6] = zero(integrator.fsalfirst)
+    integrator.k[7] = zero(integrator.fsalfirst)
     nothing
 end
 
@@ -140,6 +144,23 @@ function initialize!(integrator, cache::RadauIIA9Cache)
     nothing
 end
 
+function initialize!(integrator, cache::AdaptiveRadauConstantCache)
+    max_stages = (integrator.alg.max_order - 1) ÷ 4 * 2 + 1
+    integrator.kshortsize = max_stages + 2
+    integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
+    integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
+    OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
+
+    # Avoid undefined entries if k is an array of arrays
+    integrator.fsallast = zero(integrator.fsalfirst)
+    integrator.k[1] = integrator.fsalfirst
+    integrator.k[2] = integrator.fsallast
+    for i in 3 : max_stages + 2
+        integrator.k[i] = zero(integrator.fsallast)
+    end
+    nothing
+end
+
 function initialize!(integrator, cache::AdaptiveRadauCache)
     max_stages = (integrator.alg.max_order - 1) ÷ 4 * 2 + 1
     integrator.kshortsize = max_stages + 2
@@ -616,8 +637,8 @@ end
         if alg.extrapolant != :constant
             integrator.k[3] = (z2 - z3)/c2m1
             tmp = @.. (z1 - z2)/c1mc2
-            integrator.k[4] = (tmp - cache.cont1)/c1m1
-            integrator.k[5] = cache.cont2-(tmp - z1 / c1) / c2
+            integrator.k[4] = (tmp - integrator.k[3])/c1m1
+            integrator.k[5] = integrator.k[4]-(tmp - z1 / c1) / c2
         end
     end
 

lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl

@@ -9,7 +9,7 @@ for prob in [prob_ode_linear, prob_ode_2Dlinear]
     @test sim21.𝒪est[:L2]≈4 atol=testTol
 end
 
-sim21 = test_convergence(1 ./ 2 .^ (2.5:-1:0.5), prob, RadauIIA9(), dense_errors = true)
+sim21 = test_convergence(1 ./ 2 .^ (2.5:-1:0.5), prob_ode_linear, RadauIIA9(), dense_errors = true)
 @test sim21.𝒪est[:final]≈8 atol=testTol
 @test sim21.𝒪est[:L2]≈6 atol=testTol
 
@@ -32,9 +32,6 @@ for i in [5, 9, 13, 17, 21, 25], prob in [prob_ode_linear_big, prob_ode_2Dlinear
     @test sim21.𝒪est[:L2] ≈ ((i + 3) ÷ 2) atol=testTol
 end
 
-dts = 1 ./ 2 .^ (4.25:-1:0.25)
-local sim21 = test_convergence(dts, prob_ode_2Dlinear_big, AdaptiveRadau(min_order = 5, max_order = 5), dense_errors = true)
-
 #threaded tests
 using OrdinaryDiffEqCore
 for i in [5, 9, 13, 17, 21, 25], prob in [prob_ode_linear_big, prob_ode_2Dlinear_big]