Merge #130

130: Add convergence tests with ODEConvergenceTester r=charleskawczynski a=charleskawczynski

This PR adds convergence tests using ODEConvergenceTester.jl, which yields the following output:

<img width="500" alt="Screen Shot 2022-12-26 at 11 01 44 AM" src="https://user-images.githubusercontent.com/1880641/209576952-585a2ac9-7bfa-4267-9d01-d2403e26c739.png">

This was added into the test suite for now. I'd like the table to be a part of the docs, and tests in the test suite. Before doing this, we can improve the test problem documentation, since convergence orders for "`ark_analytic`" demonstrates both failed and super convergence, which suggests that there's something wrong or special about this problem compared to the others. The 903 super-convergence is especially striking.

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

test/Project.toml

@@ -14,7 +14,9 @@ Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearOperators = "5c8ed15e-5a4c-59e4-a42b-c7e8811fb125"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
+ODEConvergenceTester = "42a5c2e1-f365-4540-8ca5-3684de3ecd95"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
@@ -23,3 +25,5 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 OrdinaryDiffEq = "5.71.2"
+PrettyTables = "2"
+ODEConvergenceTester = "0.2"

test/convergence.jl

@@ -1,4 +1,5 @@
 using ClimaTimeSteppers, LinearAlgebra, Test
+import PrettyTables
 
 include(joinpath(@__DIR__, "convergence_utils.jl"))
 include(joinpath(@__DIR__, "utils.jl"))
@@ -41,6 +42,43 @@ include(joinpath(@__DIR__, "problems.jl"))
     end
 end
 
+#=
+using Revise; include("test/convergence.jl")
+results = tabulate_convergence_orders();
+=#
+function tabulate_convergence_orders()
+    tabs = [
+        ARS111,
+        ARS121,
+        ARS122,
+        ARS232,
+        ARS222,
+        IMKG232a,
+        IMKG232b,
+        IMKG242a,
+        IMKG242b,
+        IMKG252a,
+        IMKG252b,
+        IMKG253a,
+        IMKG253b,
+        IMKG254a,
+        IMKG254b,
+        IMKG254c,
+        HOMMEM1,
+        ARS233,
+        ARS343,
+        ARS443,
+        IMKG342a,
+        IMKG343a,
+        DBM453,
+    ]
+    test_cases = all_test_cases(Float64)
+    results = convergence_order_results(tabs, test_cases)
+    tabulate_convergence_orders(test_cases, tabs, results)
+    return results
+end
+tabulate_convergence_orders()
+
 #=
 if ArrayType == Array
 for (prob, sol) in [

test/convergence_utils.jl

@@ -1,3 +1,5 @@
+import ODEConvergenceTester as OCT
+import OrdinaryDiffEq as ODE
 
 """
     DirectSolver
@@ -47,3 +49,84 @@ function convergence_order(prob, sol, method, dts; kwargs...)
     _, order_est = hcat(ones(length(dts)), log2.(dts)) \ log2.(errs)
     return order_est
 end
+
+function default_expected_order(alg, tab)
+    return if tab isa CTS.AbstractIMEXARKTableau
+        CTS.theoretical_convergence_order(tab)
+    else
+        ODE.alg_order(alg)
+    end
+end
+
+function test_convergence_order!(test_case, tab, results = Dict(); refinement_range)
+    prob, alg = problem_algo(test_case, tab)
+    expected_order = default_expected_order(alg, tab())
+    cr = OCT.refinement_study(
+        prob,
+        alg;
+        verbose = false,
+        expected_order,
+        refinement_range, # ::UnitRange, 2:4 is more fine than 1:3
+    )
+    computed_order = maximum(cr.computed_order)
+    results[tab, test_case.test_name] = (; expected_order, computed_order)
+    return nothing
+end
+
+distance(theoretic, computed) = abs(computed - theoretic) / theoretic
+pass_conv(theoretic, computed) = distance(theoretic, computed) * 100 < 10
+fail_conv(theoretic, computed) = !pass_conv(computed, theoretic) && !super_conv(theoretic, computed)
+super_conv(theoretic, computed) = (computed - theoretic) / theoretic * 100 > 10
+
+#= Calls `test_convergence_order!` for each combination of test case
+and algorithm, returns a `Dict` of the results. =#
+function convergence_order_results(tabs, test_cases)
+    results = Dict()
+    for test_case in test_cases
+        @info "------------------ Test case: $(test_case.test_name)"
+        for tab in tabs
+            # @info "Running refinement study on $(nameof(tab))"
+            test_convergence_order!(test_case, tab, results; refinement_range = 5:9)
+        end
+    end
+    return results
+end
+
+function tabulate_convergence_orders(test_cases, tabs, results)
+    columns = map(test_cases) do test_case
+        map(tab -> results[tab, test_case.test_name], tabs)
+    end
+    expected_order = map(tab -> default_expected_order(nothing, tab()), tabs)
+    tab_names = map(tab -> "$tab ($(default_expected_order(nothing, tab())))", tabs)
+    data = hcat(columns...)
+    summary(result) = result.computed_order
+    data_summary = map(d -> summary(d), data)
+
+    table_data = hcat(tab_names, data_summary)
+    precentage_fail = sum(fail_conv.(getindex.(data, 1), getindex.(data, 2))) / length(data) * 100
+    @info "Percentage of failed convergence order tests: $precentage_fail"
+    fail_conv_hl = PrettyTables.Highlighter(
+        (data, i, j) -> j ≠ 1 && fail_conv(expected_order[i], data[i, j]),
+        PrettyTables.crayon"red bold",
+    )
+    super_conv_hl = PrettyTables.Highlighter(
+        (data, i, j) -> j ≠ 1 && super_conv(expected_order[i], data[i, j]),
+        PrettyTables.crayon"yellow bold",
+    )
+    tab_column_hl = PrettyTables.Highlighter((data, i, j) -> j == 1, PrettyTables.crayon"green bold")
+    test_case_names = map(test_case -> test_case.test_name, test_cases)
+
+    header = (["Tableau (theoretic)", test_case_names...],
+    # ["", ["" for tc in test_case_names]...],
+    )
+
+    PrettyTables.pretty_table(
+        table_data;
+        header_crayon = PrettyTables.crayon"green bold",
+        highlighters = (tab_column_hl, fail_conv_hl, super_conv_hl),
+        title = "Computed convergence orders, red=fail, yellow=super-convergence",
+        header,
+        alignment = :c,
+        crop = :none,
+    )
+end

test/problems.jl

@@ -428,3 +428,14 @@ function ark_analytic_sys_test_cts(::Type{FT}) where {FT}
         Wfact! = (W, Y, _, Δt, t) -> W .= Δt .* A .- I,
     )
 end
+
+function all_test_cases(::Type{FT}) where {FT}
+    return [
+        ark_analytic_nonlin_test_cts(FT),
+        ark_analytic_sys_test_cts(FT),
+        ark_analytic_test_cts(FT),
+        # onewaycouple_mri_test_cts(FT),
+        # climacore_2Dheat_test_cts(FT),
+        # climacore_1Dheat_test_cts(FT),
+    ]
+end