Formatted with Runic

Author: Sreeram-Shankar

lib/OrdinaryDiffEqFIRK/src/alg_utils.jl

@@ -33,4 +33,4 @@ alg_adaptive_order(alg::GaussLegendre) = 2 * alg.num_stages
 has_stiff_interpolation(::GaussLegendre) = true
 
 get_current_alg_order(alg::GaussLegendre, cache) = 2 * alg.num_stages
-get_current_adaptive_order(alg::GaussLegendre, cache) = 2 * alg.num_stages
+get_current_adaptive_order(alg::GaussLegendre, cache) = 2 * alg.num_stages

lib/OrdinaryDiffEqFIRK/src/algorithms.jl

@@ -302,4 +302,4 @@ function GaussLegendre(;
         num_stages,
         AD_choice
     )
-end
+end

lib/OrdinaryDiffEqFIRK/src/firk_caches.jl

@@ -861,4 +861,4 @@ function alg_cache(
         tmp, atmp, jac_config, linsolve, rtol, atol, dt, dt,
         Convergence, alg.step_limiter!, num_stages
     )
-end
+end

lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl

@@ -2257,9 +2257,13 @@ function initialize!(integrator, cache::GaussLegendreConstantCache)
 
     # adaptive Richardson controller requires num_stages >= 2
     if integrator.opts.adaptive && integrator.alg.num_stages < 2
-        throw(ArgumentError("GaussLegendre with num_stages = $(integrator.alg.num_stages) " *
-                            "does not support adaptive stepping (Richardson controller " *
-                            "requires num_stages ≥ 2). Use num_stages ≥ 2 or pass adaptive = false."))
+        throw(
+            ArgumentError(
+                "GaussLegendre with num_stages = $(integrator.alg.num_stages) " *
+                    "does not support adaptive stepping (Richardson controller " *
+                    "requires num_stages ≥ 2). Use num_stages ≥ 2 or pass adaptive = false."
+            )
+        )
     end
     return nothing
 end
@@ -2276,9 +2280,13 @@ function initialize!(integrator, cache::GaussLegendreCache)
     OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
     if integrator.opts.adaptive
         if cache.num_stages < 2
-            throw(ArgumentError("GaussLegendre with num_stages = $(cache.num_stages) " *
-                                "does not support adaptive stepping (Richardson controller " *
-                                "requires num_stages ≥ 2). Use num_stages ≥ 2 or pass adaptive = false."))
+            throw(
+                ArgumentError(
+                    "GaussLegendre with num_stages = $(cache.num_stages) " *
+                        "does not support adaptive stepping (Richardson controller " *
+                        "requires num_stages ≥ 2). Use num_stages ≥ 2 or pass adaptive = false."
+                )
+            )
         end
         (; abstol, reltol) = integrator.opts
         if reltol isa Number
@@ -2346,10 +2354,13 @@ end
         end
 
         ndwprev = ndw
-        ndw = sum(internalnorm(
-                      calculate_residuals(dw[i], uprev_local, uprev_local, atol, rtol, internalnorm, t_local),
-                      t_local)
-                  for i in 1:num_stages)
+        ndw = sum(
+            internalnorm(
+                    calculate_residuals(dw[i], uprev_local, uprev_local, atol, rtol, internalnorm, t_local),
+                    t_local
+                )
+                for i in 1:num_stages
+        )
 
         if iter > 1
             θ = ndw / ndwprev
@@ -2390,8 +2401,10 @@ end
     return (u_out, z, true)
 end
 
-@muladd function perform_step!(integrator, cache::GaussLegendreConstantCache,
-        repeat_step = false)
+@muladd function perform_step!(
+        integrator, cache::GaussLegendreConstantCache,
+        repeat_step = false
+    )
     (; t, dt, uprev, f, p) = integrator
     (; num_stages) = cache
     (; internalnorm, abstol, reltol, adaptive) = integrator.opts
@@ -2405,7 +2418,8 @@ end
     if adaptive && num_stages >= 2
         # Richardson step-doubling: one full step at dt, two successive half-steps.
         u_H, z_full, ok1 = _gausslegendre_substep_constant(
-            integrator, cache, alg, uprev, t, dt, J, atol, rtol)
+            integrator, cache, alg, uprev, t, dt, J, atol, rtol
+        )
         if !ok1
             integrator.force_stepfail = true
             integrator.stats.nnonlinconvfail += 1
@@ -2414,15 +2428,17 @@ end
 
         half_dt = dt / 2
         u_h1, _, ok2 = _gausslegendre_substep_constant(
-            integrator, cache, alg, uprev, t, half_dt, J, atol, rtol)
+            integrator, cache, alg, uprev, t, half_dt, J, atol, rtol
+        )
         if !ok2
             integrator.force_stepfail = true
             integrator.stats.nnonlinconvfail += 1
             return
         end
 
         u_h, _, ok3 = _gausslegendre_substep_constant(
-            integrator, cache, alg, u_h1, t + half_dt, half_dt, J, atol, rtol)
+            integrator, cache, alg, u_h1, t + half_dt, half_dt, J, atol, rtol
+        )
         if !ok3
             integrator.force_stepfail = true
             integrator.stats.nnonlinconvfail += 1
@@ -2432,13 +2448,15 @@ end
         p_order = 2 * num_stages
         utilde = @.. (u_h - u_H) / (2^p_order - 1)
         integrator.EEst = internalnorm(
-            calculate_residuals(utilde, uprev, u_h, atol, rtol, internalnorm, t), t)
+            calculate_residuals(utilde, uprev, u_h, atol, rtol, internalnorm, t), t
+        )
 
         u = u_h
         z_for_fsal = z_full
     else
         u_out, z_out, ok = _gausslegendre_substep_constant(
-            integrator, cache, alg, uprev, t, dt, J, atol, rtol)
+            integrator, cache, alg, uprev, t, dt, J, atol, rtol
+        )
         if !ok
             integrator.force_stepfail = true
             integrator.stats.nnonlinconvfail += 1
@@ -2471,8 +2489,10 @@ end
         u_dest, uprev_local, t_local, dt_local, J_local,
         cache::GaussLegendreCache, integrator, alg
     )
-    (; tab, κ, z, dw, ks, W, ubuff, tmp, atmp, linsolve,
-        rtol, atol, num_stages) = cache
+    (;
+        tab, κ, z, dw, ks, W, ubuff, tmp, atmp, linsolve,
+        rtol, atol, num_stages,
+    ) = cache
     (; A, b, c) = tab
     (; internalnorm) = integrator.opts
     (; maxiters) = alg
@@ -2510,8 +2530,10 @@ end
         end
 
         needfactor = iter == 1
-        linres = dolinsolve(integrator, linsolve;
-            A = needfactor ? W : nothing, b = ubuff, linu = ubuff)
+        linres = dolinsolve(
+            integrator, linsolve;
+            A = needfactor ? W : nothing, b = ubuff, linu = ubuff
+        )
         cache.linsolve = linres.cache
         integrator.stats.nsolve += 1
 
@@ -2520,10 +2542,12 @@ end
         end
 
         ndwprev = ndw
-        ndw = sum(begin
-                      calculate_residuals!(atmp, dw[i], uprev_local, uprev_local, atol, rtol, internalnorm, t_local)
-                      internalnorm(atmp, t_local)
-                  end for i in 1:num_stages)
+        ndw = sum(
+            begin
+                    calculate_residuals!(atmp, dw[i], uprev_local, uprev_local, atol, rtol, internalnorm, t_local)
+                    internalnorm(atmp, t_local)
+                end for i in 1:num_stages
+        )
 
         if iter > 1
             θ = ndw / ndwprev
@@ -2567,16 +2591,18 @@ end
 
 @muladd function perform_step!(integrator, cache::GaussLegendreCache, repeat_step = false)
     (; t, dt, uprev, u, f, p, fsallast) = integrator
-    (; atmp, J, z, z_last, u_full, u_half, rtol, atol,
-        step_limiter!, num_stages) = cache
+    (;
+        atmp, J, z, z_last, u_full, u_half, rtol, atol,
+        step_limiter!, num_stages,
+    ) = cache
     (; internalnorm, adaptive) = integrator.opts
     alg = unwrap_alg(integrator, true)
 
     new_jac = do_newJ(integrator, alg, cache, repeat_step)
     new_jac && (calc_J!(J, integrator, cache); cache.W_γdt = dt)
 
     if adaptive && num_stages >= 2
-        # fll step at dt 
+        # fll step at dt
         ok1 = _gausslegendre_substep!(u_full, uprev, t, dt, J, cache, integrator, alg)
         if !ok1
             integrator.force_stepfail = true
@@ -2609,7 +2635,7 @@ end
 
         p_order = 2 * num_stages
         denom = 2^p_order - 1
-        @.. u_full = (u - u_full) / denom       
+        @.. u_full = (u - u_full) / denom
         calculate_residuals!(atmp, u_full, uprev, u, atol, rtol, internalnorm, t)
         integrator.EEst = internalnorm(atmp, t)
     else
@@ -2637,4 +2663,4 @@ end
     f(fsallast, u, p, t + dt)
     OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
     return
-end
+end

lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl

@@ -427,4 +427,4 @@ const GaussLegendreTableauCache = Dict{
     (Float64, Float64, 3) => generateGaussLegendreTableau(Float64, Float64, 3),
     (Float64, Float64, 4) => generateGaussLegendreTableau(Float64, Float64, 4),
     (Float64, Float64, 5) => generateGaussLegendreTableau(Float64, Float64, 5),
-)
+)

lib/OrdinaryDiffEqFIRK/test/gausslegendre_adaptive_tests.jl

@@ -1,16 +1,18 @@
 using OrdinaryDiffEqFIRK, DiffEqDevTools, Test, LinearAlgebra
 import ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear
 
-testTol = 0.6  
+testTol = 0.6
 
 # test orders 2 and 3 for convergence with fixed dt, anything higher is too sensitive to floating point precision
 @testset "GaussLegendre: fixed-dt empirical order (s = 2, 3)" begin
     dts = Float64.(1 ./ 2 .^ (5:-1:2))  # 1/32 … 1/4
     for s in 2:3
         alg = GaussLegendre(num_stages = s; maxiters = 100)
-        sim = test_convergence(dts, prob_ode_linear, alg;
+        sim = test_convergence(
+            dts, prob_ode_linear, alg;
             dense_errors = false,
-            abstol = 1e-12, reltol = 1e-12)
+            abstol = 1.0e-12, reltol = 1.0e-12
+        )
         @test sim.𝒪est[:final] ≈ 2 * s atol = testTol
     end
 end
@@ -19,45 +21,57 @@ end
 @testset "GaussLegendre: fixed-dt accuracy (s = 4, order 8)" begin
     s = 4
     alg = GaussLegendre(num_stages = s; maxiters = 100)
-    sol = solve(prob_ode_linear, alg; adaptive = false, dt = 1 // 256,
-        abstol = 1e-14, reltol = 1e-14)
+    sol = solve(
+        prob_ode_linear, alg; adaptive = false, dt = 1 // 256,
+        abstol = 1.0e-14, reltol = 1.0e-14
+    )
     @test sol.retcode == ReturnCode.Success
     exact = prob_ode_linear.u0 * exp(1.01 * (sol.t[end] - sol.t[1]))
-    @test isapprox(sol.u[end], exact; rtol = 1e-9, atol = 1e-12)
+    @test isapprox(sol.u[end], exact; rtol = 1.0e-9, atol = 1.0e-12)
 end
 
 # test adaptive stepping with tolerance matching
 @testset "GaussLegendre: adaptive run matches tolerance" begin
     for s in 2:4
-        reltol = 1e-6
-        abstol = 1e-9
-        sol = solve(prob_ode_linear, GaussLegendre(num_stages = s);
-            reltol = reltol, abstol = abstol)
+        reltol = 1.0e-6
+        abstol = 1.0e-9
+        sol = solve(
+            prob_ode_linear, GaussLegendre(num_stages = s);
+            reltol = reltol, abstol = abstol
+        )
         @test sol.retcode == ReturnCode.Success
         exact = prob_ode_linear.u0 * exp(1.01 * (sol.t[end] - sol.t[1]))
-        @test isapprox(sol.u[end], exact; rtol = 1e-3, atol = 1e-6)
+        @test isapprox(sol.u[end], exact; rtol = 1.0e-3, atol = 1.0e-6)
     end
 end
 
 # test that Richardson step-doubling tightens the step count when the tolerance is tightened
 @testset "GaussLegendre: Richardson tightens step count when tol tightens" begin
     s = 3
-    sol_loose = solve(prob_ode_linear, GaussLegendre(num_stages = s);
-        reltol = 1e-3, abstol = 1e-6)
-    sol_tight = solve(prob_ode_linear, GaussLegendre(num_stages = s);
-        reltol = 1e-8, abstol = 1e-10)
+    sol_loose = solve(
+        prob_ode_linear, GaussLegendre(num_stages = s);
+        reltol = 1.0e-3, abstol = 1.0e-6
+    )
+    sol_tight = solve(
+        prob_ode_linear, GaussLegendre(num_stages = s);
+        reltol = 1.0e-8, abstol = 1.0e-10
+    )
     @test length(sol_tight.t) >= length(sol_loose.t)
 end
 
 # test that num_stages = 1 with adaptive throws
 @testset "GaussLegendre: num_stages = 1 with adaptive throws" begin
-    @test_throws ArgumentError solve(prob_ode_linear, GaussLegendre(num_stages = 1);
-        reltol = 1e-4, abstol = 1e-7)
+    @test_throws ArgumentError solve(
+        prob_ode_linear, GaussLegendre(num_stages = 1);
+        reltol = 1.0e-4, abstol = 1.0e-7
+    )
 end
 
 # test that num_stages = 1 with adaptive = false runs
 @testset "GaussLegendre: num_stages = 1 with adaptive = false runs" begin
-    sol = solve(prob_ode_linear, GaussLegendre(num_stages = 1);
-        adaptive = false, dt = 1 // 32)
+    sol = solve(
+        prob_ode_linear, GaussLegendre(num_stages = 1);
+        adaptive = false, dt = 1 // 32
+    )
     @test sol.retcode == ReturnCode.Success
-end
+end