Merge #97

97: Modularize remaining unit tests r=charleskawczynski a=charleskawczynski

This PR:
 - Changes the `test/problems.jl` script to instead define functions
 - Qualifies some methods
 - Moves more tests into safetestsets
 - Includes the utility functions only in the portion of the unit tests where needed.

Co-authored-by: Charles Kawczynski <[email protected]>

Author: bors[bot]

docs/src/algorithms.md

@@ -25,8 +25,6 @@ The convergence orders of the provided methods are verified using test cases fro
 using Pkg
 Pkg.activate("../../test")
 Pkg.instantiate()
-include("../../test/problems.jl")
-include("../../test/utils.jl")
 include("../../test/convergence.jl")
 Pkg.activate(".")
 ```

perf/flame.jl

@@ -24,9 +24,9 @@ function do_work!(integrator, not_generated_cache)
 end
 problem_str = parsed_args["problem"]
 prob = if problem_str=="ode_fun"
-    split_linear_prob_wfact_split
+    split_linear_prob_wfact_split()
 elseif problem_str=="fe"
-    split_linear_prob_wfact_split_fe
+    split_linear_prob_wfact_split_fe()
 else
     error("Bad option")
 end

perf/jet.jl

@@ -24,9 +24,9 @@ function config_integrators(problem)
     return (; integrator_generated=integrator, not_generated_integrator)
 end
 prob = if parsed_args["problem"]=="ode_fun"
-    split_linear_prob_wfact_split
+    split_linear_prob_wfact_split()
 elseif parsed_args["problem"]=="fe"
-    split_linear_prob_wfact_split_fe
+    split_linear_prob_wfact_split_fe()
 else
     error("Bad option")
 end

test/compare_generated.jl

@@ -6,8 +6,8 @@ include("problems.jl")
     algorithm = ARS343(NewtonsMethod(; max_iters = 2))
     dt = 0.01
     for problem in (
-        split_linear_prob_wfact_split,
-        split_linear_prob_wfact_split_fe,
+        split_linear_prob_wfact_split(),
+        split_linear_prob_wfact_split_fe(),
     )
         integrator = DiffEqBase.init(deepcopy(problem), algorithm; dt)
         not_generated_integrator = deepcopy(integrator)

test/convergence.jl

@@ -1,40 +1,46 @@
 using ClimaTimeSteppers, LinearAlgebra, Test
 
-dts = 0.5 .^ (4:7)
-
-for (prob, sol, tscale) in [
-    (linear_prob, linear_sol, 1)
-    (sincos_prob, sincos_sol, 1)
-]
-
-    @test convergence_order(prob, sol, LSRKEulerMethod(), dts.*tscale)     ≈ 1 atol=0.1
-    @test convergence_order(prob, sol, LSRK54CarpenterKennedy(), dts.*tscale)       ≈ 4 atol=0.05
-    @test convergence_order(prob, sol, LSRK144NiegemannDiehlBusch(),  dts.*tscale)   ≈ 4 atol=0.05
+include(joinpath(@__DIR__, "convergence_utils.jl"))
+include(joinpath(@__DIR__, "utils.jl"))
+include(joinpath(@__DIR__, "problems.jl"))
+
+@testset "LSRK and SSP convergence" begin
+    dts = 0.5 .^ (4:7)
+
+    for (prob, sol, tscale) in [
+        (linear_prob(), linear_sol, 1)
+        (sincos_prob(), sincos_sol, 1)
+    ]
+
+        @test convergence_order(prob, sol, LSRKEulerMethod(), dts.*tscale)     ≈ 1 atol=0.1
+        @test convergence_order(prob, sol, LSRK54CarpenterKennedy(), dts.*tscale)       ≈ 4 atol=0.05
+        @test convergence_order(prob, sol, LSRK144NiegemannDiehlBusch(),  dts.*tscale)   ≈ 4 atol=0.05
+
+        @test convergence_order(prob, sol, SSPRK22Heuns(), dts.*tscale)                 ≈ 2 atol=0.05
+        @test convergence_order(prob, sol, SSPRK22Ralstons(), dts.*tscale)              ≈ 2 atol=0.05
+        @test convergence_order(prob, sol, SSPRK33ShuOsher(), dts.*tscale)              ≈ 3 atol=0.05
+        @test convergence_order(prob, sol, SSPRK34SpiteriRuuth(), dts.*tscale)          ≈ 3 atol=0.05
+    end
 
-    @test convergence_order(prob, sol, SSPRK22Heuns(), dts.*tscale)                 ≈ 2 atol=0.05
-    @test convergence_order(prob, sol, SSPRK22Ralstons(), dts.*tscale)              ≈ 2 atol=0.05
-    @test convergence_order(prob, sol, SSPRK33ShuOsher(), dts.*tscale)              ≈ 3 atol=0.05
-    @test convergence_order(prob, sol, SSPRK34SpiteriRuuth(), dts.*tscale)          ≈ 3 atol=0.05
-end
+    for (prob, sol, tscale) in [
+        (linear_prob_wfactt(), linear_sol, 1)
+    ]
+        @test convergence_order(prob, sol, SSPKnoth(linsolve=linsolve_direct), dts.*tscale) ≈ 2 atol=0.05
 
-for (prob, sol, tscale) in [
-    (linear_prob_wfactt, linear_sol, 1)
-]
-    @test convergence_order(prob, sol, SSPKnoth(linsolve=linsolve_direct), dts.*tscale) ≈ 2 atol=0.05
+    end
 
-end
 
+    # ForwardEulerODEFunction
+    for (prob, sol, tscale) in [
+        (linear_prob_fe(), linear_sol, 1)
+        (sincos_prob_fe(), sincos_sol, 1)
+    ]
+        @test convergence_order(prob, sol, SSPRK22Heuns(), dts.*tscale)                 ≈ 2 atol=0.05
+        @test convergence_order(prob, sol, SSPRK22Ralstons(), dts.*tscale)              ≈ 2 atol=0.05
+        @test convergence_order(prob, sol, SSPRK33ShuOsher(), dts.*tscale)              ≈ 3 atol=0.05
+        @test convergence_order(prob, sol, SSPRK34SpiteriRuuth(), dts.*tscale)          ≈ 3 atol=0.05
 
-# ForwardEulerODEFunction
-for (prob, sol, tscale) in [
-    (linear_prob_fe, linear_sol, 1)
-    (sincos_prob_fe, sincos_sol, 1)
-]
-    @test convergence_order(prob, sol, SSPRK22Heuns(), dts.*tscale)                 ≈ 2 atol=0.05
-    @test convergence_order(prob, sol, SSPRK22Ralstons(), dts.*tscale)              ≈ 2 atol=0.05
-    @test convergence_order(prob, sol, SSPRK33ShuOsher(), dts.*tscale)              ≈ 3 atol=0.05
-    @test convergence_order(prob, sol, SSPRK34SpiteriRuuth(), dts.*tscale)          ≈ 3 atol=0.05
-
+    end
 end
 
 ENV["GKSwstype"] = "nul" # avoid displaying plots
@@ -46,14 +52,14 @@ ENV["GKSwstype"] = "nul" # avoid displaying plots
     algs2 = (algs2..., IMKG254b, IMKG254c, HOMMEM1)
     algs3 = (ARS233, ARS343, ARS443, IMKG342a, IMKG343a, DBM453)
     dict = Dict(((algs1 .=> 1)..., (algs2 .=> 2)..., (algs3 .=> 3)...))
-    test_algs("IMEX ARK", dict, ark_analytic_nonlin_test, 400)
-    test_algs("IMEX ARK", dict, ark_analytic_sys_test, 60)
+    test_algs("IMEX ARK", dict, ark_analytic_nonlin_test(Float64), 400)
+    test_algs("IMEX ARK", dict, ark_analytic_sys_test(Float64), 60)
 
     # For some bizarre reason, ARS121 converges with order 2 for ark_analytic,
     # even though it is only a 1st order method.
     dict′ = copy(dict)
     dict′[ARS121] = 2
-    test_algs("IMEX ARK", dict′, ark_analytic_test, 16000)
+    test_algs("IMEX ARK", dict′, ark_analytic_test(Float64), 16000)
 end
 
 #=

test/convergence_utils.jl

@@ -0,0 +1,49 @@
+
+"""
+    DirectSolver
+
+A linear solver which forms the full matrix of a linear operator and its LU factorization.
+"""
+struct DirectSolver end
+
+DirectSolver(args...) = DirectSolver()
+
+function (::DirectSolver)(x,A,b,matrix_updated; kwargs...)
+  n = length(x)
+  M = mapslices(y -> mul!(similar(y), A, y), Matrix{eltype(x)}(I,n,n), dims=1)
+  x .= M \ b
+end
+
+"""
+    convergence_rates(problem, solution, method, dts; kwargs...)
+
+Compute the errors rates of `method` on `problem` by comparing to `solution`
+on the set of `dt` values in `dts`. Extra `kwargs` are passed to `solve`
+
+`solution` should be a function with a method `solution(u0, p, t)`.
+"""
+function convergence_errors(prob, sol, method, dts; kwargs...)
+  errs = map(dts) do dt
+       # copy the problem so we don't mutate u0
+      prob_copy = deepcopy(prob)
+      u = solve(prob_copy, method; dt=dt, saveat=(prob.tspan[2],), kwargs...)
+      norm(u .- sol(prob.u0, prob.p, prob.tspan[end]))
+  end
+  return errs
+end
+
+
+"""
+  convergence_order(problem, solution, method, dts; kwargs...)
+
+Estimate the order of the rate of convergence of `method` on `problem` by comparing to
+`solution` the set of `dt` values in `dts`.
+
+`solution` should be a function with a method `solution(u0, p, t)`.
+"""
+function convergence_order(prob, sol, method, dts; kwargs...)
+  errs = convergence_errors(prob, sol, method, dts; kwargs...)
+  # find slope coefficient in log scale
+  _,order_est = hcat(ones(length(dts)), log2.(dts)) \ log2.(errs)
+  return order_est
+end

test/integrator.jl

@@ -4,7 +4,7 @@ import OrdinaryDiffEq
 include("problems.jl")
 
 @testset "integrator save times" begin
-    test_case = constant_tendency_test
+    test_case = constant_tendency_test(Float64)
     (; prob, analytic_sol) = test_case
     for alg in (SSPRK33ShuOsher(), OrdinaryDiffEq.SSPRK33()),
         reverse_prob in (false, true),
@@ -121,7 +121,7 @@ end
     # OrdinaryDiffEq does not save at t0′ after reinit! unless erase_sol is
     # true, so this test does not include a comparison with OrdinaryDiffEq.
     alg = SSPRK33ShuOsher()
-    test_case = constant_tendency_test
+    test_case = constant_tendency_test(Float64)
     (; prob, analytic_sol) = test_case
     for reverse_prob in (false, true)
         if reverse_prob
@@ -161,7 +161,7 @@ end
 
 @testset "integrator step past end time" begin
     alg = SSPRK33ShuOsher()
-    test_case = constant_tendency_test
+    test_case = constant_tendency_test(Float64)
     (; prob, analytic_sol) = test_case
     t0, tf = prob.tspan
     dt = tf - t0