From 9ff67b0ae80d16ad7b948707516d256ff17a46f2 Mon Sep 17 00:00:00 2001 From: Dennis Yatunin Date: Fri, 31 Jan 2025 17:41:19 -0800 Subject: [PATCH] Clean up checks for convergence and limiters --- .buildkite/pipeline.yml | 10 +- docs/src/dev/compute_convergence.jl | 275 +++++--------- docs/src/dev/limiter_analysis.jl | 46 +-- docs/src/dev/limiter_summary.jl | 131 +++++-- docs/src/dev/report_gen_alg.jl | 142 +++---- docs/src/dev/summarize_convergence.jl | 441 +++++++--------------- docs/src/plotting_utils.jl | 190 ---------- ext/ClimaTimeSteppersBenchmarkToolsExt.jl | 14 +- test/problems.jl | 36 +- 9 files changed, 444 insertions(+), 841 deletions(-) delete mode 100644 docs/src/plotting_utils.jl diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 29210987..f403aab9 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -267,6 +267,12 @@ steps: - label: "Convergence: IMKG343a" command: "julia --color=yes --check-bounds=yes --project=.buildkite docs/src/dev/report_gen_alg.jl --alg IMKG343a" + - label: "Convergence: ARK437L2SA1" + command: "julia --color=yes --check-bounds=yes --project=.buildkite docs/src/dev/report_gen_alg.jl --alg ARK437L2SA1" + + - label: "Convergence: ARK548L2SA2" + command: "julia --color=yes --check-bounds=yes --project=.buildkite docs/src/dev/report_gen_alg.jl --alg ARK548L2SA2" + - label: "Convergence: SSPKnoth" command: "julia --color=yes --check-bounds=yes --project=.buildkite docs/src/dev/report_gen_alg.jl --alg SSPKnoth" @@ -276,10 +282,10 @@ steps: - label: "Summarize convergence" command: "julia --color=yes --check-bounds=yes --project=.buildkite docs/src/dev/summarize_convergence.jl" - artifact_paths: "output/*" + artifact_paths: "output/convergence_*.png" depends_on: alg_convergence - label: "Summarize limiter analysis" command: "julia --color=yes --project=.buildkite docs/src/dev/limiter_summary.jl" - artifact_paths: "output/*" + artifact_paths: "output/limiter_*.png" depends_on: limiters_analysis diff --git a/docs/src/dev/compute_convergence.jl b/docs/src/dev/compute_convergence.jl index 20315041..1a03bea1 100644 --- a/docs/src/dev/compute_convergence.jl +++ b/docs/src/dev/compute_convergence.jl @@ -1,27 +1,9 @@ -import JLD2 -import Plots +using ClimaTimeSteppers +using InteractiveUtils: subtypes using Distributions: quantile, TDist -using Printf: @sprintf -using LaTeXStrings: latexstring -import DiffEqCallbacks -import ClimaTimeSteppers as CTS +using LinearAlgebra: norm -function get_algorithm_names() - all_subtypes(::Type{T}) where {T} = isabstracttype(T) ? vcat(all_subtypes.(subtypes(T))...) : [T] - algorithm_names = map(T -> T(), all_subtypes(ClimaTimeSteppers.AbstractAlgorithmName)) - return filter(name -> !(name isa ARK437L2SA1 || name isa ARK548L2SA2), algorithm_names) # too high order -end - -function get_imex_ssprk_algorithm_names() - all_subtypes(::Type{T}) where {T} = isabstracttype(T) ? vcat(all_subtypes.(subtypes(T))...) : [T] - algorithm_names = map(T -> T(), all_subtypes(ClimaTimeSteppers.IMEXSSPRKAlgorithmName)) - return algorithm_names -end - -function make_saving_callback(cb, u, t, integrator) - savevalType = typeof(cb(u, t, integrator)) - return DiffEqCallbacks.SavingCallback(cb, DiffEqCallbacks.SavedValues(typeof(t), savevalType)) -end +all_subtypes(::Type{T}) where {T} = isabstracttype(T) ? vcat(all_subtypes.(subtypes(T))...) : [T] """ imex_convergence_orders(algorithm_name) @@ -60,8 +42,7 @@ imex_convergence_orders(::ARK548L2SA2) = (5, 5, 5) imex_convergence_orders(::SSP22Heuns) = (2, 2, 2) imex_convergence_orders(::SSP33ShuOsher) = (3, 3, 3) imex_convergence_orders(::RK4) = (4, 4, 4) -# SSPKnoth is not really an IMEX method -imex_convergence_orders(::SSPKnoth) = (2, 2, 2) +imex_convergence_orders(::SSPKnoth) = (2, 3, 2) # Compute a confidence interval for the convergence order, returning the # estimated convergence order and its uncertainty. @@ -103,192 +84,136 @@ function (assuming that the algorithm converges). function predicted_convergence_order(algorithm_name::AbstractAlgorithmName, ode_function::AbstractClimaODEFunction) (imp_order, exp_order, combined_order) = imex_convergence_orders(algorithm_name) has_imp = !isnothing(ode_function.T_imp!) - has_exp = CTS.has_T_exp(ode_function) + has_exp = ClimaTimeSteppers.has_T_exp(ode_function) has_imp && !has_exp && return imp_order !has_imp && has_exp && return exp_order has_imp && has_exp && return combined_order return 0 end -function export_convergence_results(alg_name, test_problem, num_steps; kwargs...) - out_dict = Dict() - (; test_name) = test_problem - out_dict[string(test_name)] = Dict() - out_dict[string(test_name)][string(alg_name)] = Dict() - out_dict[string(test_name)]["args"] = (alg_name, test_problem, num_steps) - out_dict[string(test_name)]["kwargs"] = kwargs - compute_convergence!(out_dict, alg_name, test_problem, num_steps; kwargs...) - JLD2.save_object("convergence_$(alg_name)_$(test_problem.test_name).jld2", out_dict) -end - +""" + algorithm(algorithm_name, [linear_implicit]) -function compute_convergence!( - out_dict, - alg_name, +Generates an appropriate `DistributedODEAlgorithm` from an `AbstractAlgorithmName`. +For `IMEXAlgorithmNames`, `linear_implicit` must also be specified. One Newton +iteration is used for linear implicit problems, and two iterations are used for +nonlinear implicit problems. +""" +algorithm(algorithm_name, _) = algorithm(algorithm_name) +algorithm(algorithm_name::ClimaTimeSteppers.ERKAlgorithmName) = ExplicitAlgorithm(algorithm_name) +algorithm(algorithm_name::ClimaTimeSteppers.SSPKnoth) = + ClimaTimeSteppers.RosenbrockAlgorithm(ClimaTimeSteppers.tableau(ClimaTimeSteppers.SSPKnoth())) +algorithm(algorithm_name::ClimaTimeSteppers.IMEXARKAlgorithmName, linear_implicit) = + IMEXAlgorithm(algorithm_name, NewtonsMethod(; max_iters = linear_implicit ? 1 : 2)) + +rms(array) = norm(array) / sqrt(length(array)) +rms_error(u, t, sol) = rms(abs.(u .- sol(t))) + +function test_convergence!( + convergence_results, + algorithm_name, test_case, - num_steps; - num_steps_scaling_factor = 10, - order_confidence_percent = 99, - super_convergence = (), + default_num_steps; + super_convergence_algorithm_names = (), + num_steps_scaling_factor = 4, + high_order_sample_shifts = 1, numerical_reference_algorithm_name = nothing, - numerical_reference_num_steps = num_steps_scaling_factor^3 * num_steps, - full_history_algorithm_name = nothing, - average_function = array -> norm(array) / sqrt(length(array)), - average_function_str = "RMS", - only_endpoints = false, + numerical_reference_num_steps = num_steps_scaling_factor^3 * default_num_steps, + broken_tests = (), + error_on_failure = true, verbose = false, ) (; test_name, t_end, linear_implicit, analytic_sol) = test_case prob = test_case.split_prob - FT = typeof(t_end) - default_dt = t_end / num_steps - key1 = string(test_name) - key2 = string(alg_name) - - algorithm(algorithm_name::ClimaTimeSteppers.ERKAlgorithmName) = ExplicitAlgorithm(algorithm_name) - algorithm(algorithm_name::ClimaTimeSteppers.SSPKnoth) = - ClimaTimeSteppers.RosenbrockAlgorithm(ClimaTimeSteppers.tableau(ClimaTimeSteppers.SSPKnoth())) - algorithm(algorithm_name::ClimaTimeSteppers.IMEXARKAlgorithmName) = - IMEXAlgorithm(algorithm_name, NewtonsMethod(; max_iters = linear_implicit ? 1 : 2)) + default_dt = t_end / default_num_steps ref_sol = if isnothing(numerical_reference_algorithm_name) analytic_sol else - ref_alg = algorithm(numerical_reference_algorithm_name) + # TODO: Do not regenerate the reference solution for every algorithm!! ref_alg_str = string(nameof(typeof(numerical_reference_algorithm_name))) + ref_alg = algorithm(numerical_reference_algorithm_name, linear_implicit) ref_dt = t_end / numerical_reference_num_steps - verbose && - @info "Generating numerical reference solution for $test_name with $ref_alg_str (dt = $ref_dt)..." - sol = solve(deepcopy(prob), ref_alg; dt = ref_dt, save_everystep = !only_endpoints) - out_dict[string(test_name)]["numerical_ref_sol"] = sol + verbose && @info "Generating reference solution for $test_name with $ref_alg_str and dt = $ref_dt" + solve(deepcopy(prob), ref_alg; dt = ref_dt, save_everystep = true) end - - cur_avg_err(u, t) = average_function(abs.(u .- ref_sol(t))) - cur_avg_sol_and_err(u, t) = (average_function(u), average_function(abs.(u .- ref_sol(t)))) - - float_str(x) = @sprintf "%.4f" x - pow_str(x) = "10^{$(@sprintf "%.1f" log10(x))}" - function si_str(x) - if isnan(x) || x in (0, Inf, -Inf) - return string(x) - end - exponent = floor(Int, log10(x)) - mantissa = x / 10.0^exponent - return "$(float_str(mantissa)) \\times 10^{$exponent}" + numerical_reference_info = if isnothing(numerical_reference_algorithm_name) + nothing + else + ref_average_rms_error = rms(rms_error.(ref_sol.u, ref_sol.t, (analytic_sol,))) + (ref_alg_str, ref_dt, ref_average_rms_error) end - net_avg_sol_str = "\\textrm{$average_function_str}\\_\\textrm{solution}" - net_avg_err_str = "\\textrm{$average_function_str}\\_\\textrm{error}" - cur_avg_sol_str = "\\textrm{current}\\_$net_avg_sol_str" - cur_avg_err_str = "\\textrm{current}\\_$net_avg_err_str" - - linestyles = (:solid, :dash, :dot, :dashdot, :dashdotdot) - marker_kwargs = (; markershape = :circle, markeralpha = 0.5, markerstrokewidth = 0) - plot_kwargs = (; - legendposition = :outerright, - legendtitlefontpointsize = 8, - palette = :glasbey_bw_minc_20_maxl_70_n256, - size = (1000, 2000), # size in px - leftmargin = 60Plots.px, - rightmargin = 0Plots.px, - topmargin = 0Plots.px, - bottommargin = 30Plots.px, - ) - - plot1_dts = t_end ./ round.(Int, num_steps .* num_steps_scaling_factor .^ (-1:0.5:1)) - plot1 = Plots.plot(; - title = "Convergence Orders", - xaxis = (latexstring("dt"), :log10), - yaxis = (latexstring(net_avg_err_str), :log10), - legendtitle = "Convergence Order ($order_confidence_percent% CI)", - plot_kwargs..., - ) - - plot2b_min = typemax(FT) - plot2b_max = typemin(FT) - plot2a = Plots.plot(; - title = latexstring("Solutions with \$dt = $(pow_str(default_dt))\$"), - xaxis = (latexstring("t"),), - yaxis = (latexstring(cur_avg_sol_str),), - legendtitle = latexstring(net_avg_sol_str), - plot_kwargs..., - ) - plot2b = Plots.plot(; - title = latexstring("Errors with \$dt = $(pow_str(default_dt))\$"), - xaxis = (latexstring("t"),), - yaxis = (latexstring(cur_avg_err_str), :log10), - legendtitle = latexstring(net_avg_err_str), - plot_kwargs..., - ) - - cur_avg_errs_dict = Dict() - # for algorithm_name in algorithm_names - algorithm_name = alg_name - alg = algorithm(algorithm_name) alg_str = string(nameof(typeof(algorithm_name))) - predicted_order = predicted_convergence_order(algorithm_name, prob.f) - linestyle = linestyles[(predicted_order - 1) % length(linestyles) + 1] + alg = algorithm(algorithm_name, linear_implicit) + verbose && @info "Testing convergence of $alg_str for $test_name" - verbose && @info "Running $test_name with $alg_str..." - @info "Using plot1_dts=$plot1_dts" - plot1_net_avg_errs = map(plot1_dts) do plot1_dt - plot1_sol = solve(deepcopy(prob), alg; dt = plot1_dt, save_everystep = !only_endpoints) - (; u, t) = plot1_sol - cur_avg_errs = cur_avg_err.(u, t) - cur_avg_errs_dict[plot1_dt] = cur_avg_errs - verbose && @info "RMS_error(dt = $plot1_dt) = $(average_function(cur_avg_errs))" - return average_function(cur_avg_errs) + predicted_order = predicted_convergence_order(algorithm_name, prob.f) + predicted_super_convergence = algorithm_name in super_convergence_algorithm_names + num_steps_powers = (-1:0.5:1) .- high_order_sample_shifts * max(0, predicted_order - 3) / 2 + sampled_num_steps = default_num_steps .* num_steps_scaling_factor .^ num_steps_powers + sampled_dts = t_end ./ round.(Int, sampled_num_steps) + average_rms_errors = map(sampled_dts) do dt + sol = solve(deepcopy(prob), alg; dt = dt, save_everystep = true) + rms(rms_error.(sol.u, sol.t, (ref_sol,))) end - out_dict[key1][key2]["cur_avg_errs_dict"] = cur_avg_errs_dict - order, order_uncertainty = convergence_order(plot1_dts, plot1_net_avg_errs, order_confidence_percent / 100) - order_str = "$(float_str(order)) \\pm $(float_str(order_uncertainty))" - if algorithm_name in super_convergence - predicted_order += 1 - plot1_label = "$alg_str: \$$order_str\\ \\ \\ \\textbf{\\textit{SC}}\$" + verbose && @info "Sampled timesteps = $sampled_dts" + verbose && @info "Average RMS errors = $average_rms_errors" + + # Compute a 99% confidence interval for the convergence order + order, order_uncertainty = convergence_order(sampled_dts, average_rms_errors, 0.99) + verbose && @info "Convergence order = $order ± $order_uncertainty" + actual_predicted_order = predicted_order + Bool(predicted_super_convergence) + convergence_test_error = if isnan(order) + "Timestepper does not converge for $alg_str ($test_name)" + elseif abs(order - actual_predicted_order) > order_uncertainty + "Predicted order outside error bars for $alg_str ($test_name)" + elseif order_uncertainty > actual_predicted_order / 10 + "Order uncertainty too large for $alg_str ($test_name)" else - plot1_label = "$alg_str: \$$order_str\$" - end - verbose && @info "Order = $order ± $order_uncertainty" - if abs(order - predicted_order) > order_uncertainty - @warn "Predicted order outside error bars for $alg_str ($test_name)" - end - if order_uncertainty > predicted_order / 10 - @warn "Order uncertainty too large for $alg_str ($test_name)" + nothing end - - # Remove all 0s from plot2_cur_avg_errs because they cannot be plotted on a - # logarithmic scale. Record the extrema of plot2_cur_avg_errs to set ylim. - plot2_data = solve(deepcopy(prob), alg; dt = default_dt, save_everystep = true) - if any(isnan, plot2_data) - error("NaN found in plot2_data in problem $(test_name)") + if isnothing(convergence_test_error) + @assert !(algorithm_name in broken_tests) + elseif error_on_failure && !(algorithm_name in broken_tests) + error(convergence_test_error) + else + @warn convergence_test_error end - (; u, t) = plot2_data - cur_sols_and_errs = cur_avg_sol_and_err.(u, t) - out_dict[key1][key2]["plot2_data"] = (; u = cur_sols_and_errs, t) - if !isnothing(full_history_algorithm_name) - history_alg = algorithm(full_history_algorithm_name) - history_alg_name = string(nameof(typeof(full_history_algorithm_name))) - history_solve_sol = solve(deepcopy(prob), history_alg; dt = default_dt, save_everystep = true) - (; u, t) = history_solve_sol - history_solve_results = map(X -> X[1] .- ref_sol(X[2]), zip(u, t)) - history_solve_results = (; u = history_solve_results, t) - out_dict[key1][key2]["history_solve_results"] = history_solve_results - end - return out_dict + default_dt_sol = solve(deepcopy(prob), alg; dt = default_dt, save_everystep = true) + default_dt_times = default_dt_sol.t + default_dt_solutions = rms.(default_dt_sol.u) + default_dt_errors = rms_error.(default_dt_sol.u, default_dt_sol.t, (ref_sol,)) + + convergence_results[test_name] = Dict() + convergence_results[test_name]["default_dt"] = default_dt + convergence_results[test_name]["numerical_reference_info"] = numerical_reference_info + convergence_results[test_name]["all_alg_results"] = Dict() + convergence_results[test_name]["all_alg_results"][alg_str] = Dict() + alg_results = convergence_results[test_name]["all_alg_results"][alg_str] + alg_results["predicted_order"] = predicted_order + alg_results["predicted_super_convergence"] = predicted_super_convergence + alg_results["sampled_dts"] = sampled_dts + alg_results["average_rms_errors"] = average_rms_errors + alg_results["default_dt_times"] = default_dt_times + alg_results["default_dt_solutions"] = default_dt_solutions + alg_results["default_dt_errors"] = default_dt_errors + return convergence_results end -function test_unconstrained_vs_ssp_without_limiters(alg_name, test_case, num_steps) +function test_unconstrained_vs_ssp_without_limiters(algorithm_name, test_case, num_steps) prob = test_case.split_prob dt = test_case.t_end / num_steps newtons_method = NewtonsMethod(; max_iters = test_case.linear_implicit ? 1 : 2) - algorithm = IMEXAlgorithm(alg_name, newtons_method) - reference_algorithm = IMEXAlgorithm(alg_name, newtons_method, Unconstrained()) + algorithm = IMEXAlgorithm(algorithm_name, newtons_method) + reference_algorithm = IMEXAlgorithm(algorithm_name, newtons_method, Unconstrained()) solution = solve(deepcopy(prob), algorithm; dt).u[end] reference_solution = solve(deepcopy(prob), reference_algorithm; dt).u[end] - if norm(solution .- reference_solution) / norm(reference_solution) > 30 * eps(Float64) - alg_str = string(nameof(typeof(alg_name))) - @warn "Unconstrained and SSP versions of $alg_str \ - give different results for $(test_case.test_name)" + relative_error = norm(solution .- reference_solution) / norm(reference_solution) + if relative_error > 100 * eps(Float64) + error("Unconstrained and SSP versions of $algorithm_name give \ + different results for $(test_case.test_name): relative \ + error = $(round(Int, relative_error / eps(Float64))) * eps") end end diff --git a/docs/src/dev/limiter_analysis.jl b/docs/src/dev/limiter_analysis.jl index 5f0f5008..7f8e0379 100644 --- a/docs/src/dev/limiter_analysis.jl +++ b/docs/src/dev/limiter_analysis.jl @@ -1,23 +1,9 @@ -#= -julia --project=docs -using Revise; include(joinpath("docs", "src", "dev", "limiter_analysis.jl")) - -Tested algs: - - - `--alg SSP333 --use_limiter true --use_hyperdiffusion true` - - `--alg ARS343 --use_limiter true --use_hyperdiffusion true` - - `--alg SSP333 --use_limiter false --use_hyperdiffusion true` - - `--alg ARS343 --use_limiter false --use_hyperdiffusion true` - - `--alg SSP333 --use_limiter true --use_hyperdiffusion false` - - `--alg ARS343 --use_limiter true --use_hyperdiffusion false` - - `--alg SSP333 --use_limiter false --use_hyperdiffusion false` - - `--alg ARS343 --use_limiter false --use_hyperdiffusion false` -=# using ClimaTimeSteppers -using Test -using InteractiveUtils: subtypes - import ArgParse +import JLD2 + +include(joinpath(pkgdir(ClimaTimeSteppers), "test", "problems.jl")) + function parse_commandline() s = ArgParse.ArgParseSettings() ArgParse.@add_arg_table s begin @@ -37,34 +23,26 @@ function parse_commandline() parsed_args = ArgParse.parse_args(ARGS, s) return parsed_args end - -ENV["GKSwstype"] = "nul" # avoid displaying plots - -import ClimaTimeSteppers as CTS -include(joinpath(pkgdir(ClimaTimeSteppers), "test", "problems.jl")) -using ClimaTimeSteppers -import JLD2 - parsed_args = parse_commandline() -alg_name = getproperty(CTS, Symbol(parsed_args["alg"]))() +alg_str = parsed_args["alg"] use_limiter = parsed_args["use_limiter"] use_hyperdiffusion = parsed_args["use_hyperdiffusion"] - +alg_name = getproperty(ClimaTimeSteppers, Symbol(alg_str))() # This also runs with num_steps = 1000, but with larger under/overshoots; 4800 # is the value used in the paper. - -FT = Float64 num_steps = 4800 -out_dict = Dict() -test_case = horizontal_deformational_flow_test(FT; use_limiter, use_hyperdiffusion) +test_case = horizontal_deformational_flow_test(Float64; use_limiter, use_hyperdiffusion) prob = test_case.split_prob dt = test_case.t_end / num_steps algorithm = alg_name isa ClimaTimeSteppers.IMEXARKAlgorithmName ? IMEXAlgorithm(alg_name, NewtonsMethod()) : ExplicitAlgorithm(alg_name) solution = solve(deepcopy(prob), algorithm; dt).u -out_dict[alg_name, use_hyperdiffusion, use_limiter] = solution -JLD2.save_object("limiter_$(alg_name)_hyperdiff_$(use_hyperdiffusion)_lim_$(use_limiter).jld2", out_dict) +limiter_results = Dict() +limiter_results[alg_str, use_hyperdiffusion, use_limiter] = solution +mkpath("output") +filename = "limiter_$(alg_str)_hyperdiff_$(use_hyperdiffusion)_lim_$(use_limiter).jld2" +JLD2.save_object(joinpath("output", filename), limiter_results) diff --git a/docs/src/dev/limiter_summary.jl b/docs/src/dev/limiter_summary.jl index f13c1017..b54b718b 100644 --- a/docs/src/dev/limiter_summary.jl +++ b/docs/src/dev/limiter_summary.jl @@ -1,29 +1,114 @@ -#= -julia --project=docs -using Revise; include(joinpath("docs", "src", "dev", "limiter_summary.jl")) +import Plots +import JLD2 +using ClimaCorePlots +using LinearAlgebra: norm +using PrettyTables: pretty_table, ft_printf -=# -using ClimaTimeSteppers -using Test +title_str(name) = titlecase(replace(string(name), '_' => ' ')) + +# Generates the Table 1 from +# "Optimization-based limiters for the spectral element method" by Guba et al., +# and also plots the values used to generate the table. +function limiter_summary(sol_dicts, alg_strs) + table_rows = [] + mkpath("output") + for alg_str in alg_strs + plots = [] + plot_kwargs = (; + clims = (0, 1), + color = :diverging_rainbow_bgymr_45_85_c67_n256, + colorbar = false, + guide = "", + margin = 10Plots.px, + ) + for use_limiter in (false, true), use_hyperdiffusion in (false, true) + solution = sol_dicts[alg_str, use_hyperdiffusion, use_limiter] + initial_q = solution[1].ρq ./ solution[1].ρ + final_q = solution[end].ρq ./ solution[end].ρ + names = propertynames(initial_q) -import ClimaTimeSteppers as CTS -include(joinpath(pkgdir(ClimaTimeSteppers), "test", "problems.jl")) -include(joinpath(@__DIR__, "compute_convergence.jl")) -include(joinpath(@__DIR__, "..", "plotting_utils.jl")) -using ClimaTimeSteppers + if isempty(plots) + for name in names + push!(plots, Plots.plot(initial_q.:($name); plot_kwargs..., title = title_str(name))) + end + end + for name in names + push!(plots, Plots.plot(final_q.:($name); plot_kwargs..., title = "")) + end -let # Convergence - # NOTE: Some imperfections in the convergence order for SSPKnoth are to be - # expected because we are not using the exact Jacobian - vector_of_dicts = [] - files = readdir() - filter!(x -> endswith(x, ".jld2"), files) - filter!(x -> startswith(basename(x), "limiter_"), files) - for f in files - push!(vector_of_dicts, JLD2.load_object(f)) + for name in names + ϕ₀ = initial_q.:($name) + ϕ = final_q.:($name) + Δϕ₀ = maximum(ϕ₀) - minimum(ϕ₀) + ϕ_error = ϕ .- ϕ₀ + table_row = [ + alg_str;; + string(use_limiter);; + string(use_hyperdiffusion);; + title_str(name);; + max(0, -(minimum(ϕ) - minimum(ϕ₀)) / Δϕ₀);; + max(0, (maximum(ϕ) - maximum(ϕ₀)) / Δϕ₀);; + map(p -> norm(ϕ_error, p) / norm(ϕ₀, p), (1, 2, Inf))... + ] + push!(table_rows, table_row) + end + end + colorbar_plot = Plots.scatter( + [0]; + plot_kwargs..., + colorbar = true, + framestyle = :none, + legend_position = :none, + margin = 0Plots.px, + markeralpha = 0, + zcolor = [0], + ) + plot = Plots.plot( + plots..., + colorbar_plot; + layout = (Plots.@layout [Plots.grid(5, 3) a{0.1w}]), + plot_title = "Tracer specific humidity for $alg_str (Initial, \ + Final, Final w/ Hyperdiffusion, Final w/ Limiter, \ + Final w/ Hyperdiffusion & Limiter)", + size = (1600, 2000), + ) + Plots.savefig(plot, joinpath("output", "limiter_summary_$alg_str.png")) end + table = pretty_table( + vcat(table_rows...); + header = [ + "Algorithm", + "Limiter", + "Hyperdiffusion", + "Tracer Name", + "Max Undershoot", + "Max Overshoot", + "1-Norm Error", + "2-Norm Error", + "∞-Norm Error", + ], + crop = :none, + body_hlines = collect(3:3:(length(table_rows) - 1)), + formatters = ft_printf("%.4e"), + ) + println(table) +end - sol_dicts = merge(vector_of_dicts...) - algorithm_names = map(x -> x[1], collect(keys(sol_dicts))) - limiter_summary(sol_dicts, algorithm_names) +limiter_results_filenames = filter(readdir("output"; join = true)) do filename + startswith(basename(filename), "limiter_") && endswith(filename, ".jld2") end +limiter_results = merge(map(JLD2.load_object, limiter_results_filenames)...) +alg_strs = unique(map(key -> key[1], collect(keys(limiter_results)))) +limiter_summary(limiter_results, alg_strs) + +#= TO ANALYZE LIMITERS WITH ARS343 AND SSP333: +for alg_str in ("ARS343", "SSP333"), use_limiter in (false, true), use_hyperdiffusion in (false, true) + empty!(ARGS) + push!(ARGS, "--alg", alg_str) + push!(ARGS, "--use_limiter", string(use_limiter)) + push!(ARGS, "--use_hyperdiffusion", string(use_hyperdiffusion)) + @info "Running deformational flow with $alg_str, limiter = $use_limiter, hyperdif = $use_hyperdiffusion" + include("docs/src/dev/limiter_analysis.jl") +end +include("docs/src/dev/limiter_summary.jl") +=# diff --git a/docs/src/dev/report_gen_alg.jl b/docs/src/dev/report_gen_alg.jl index fbb022d8..3100cfd9 100644 --- a/docs/src/dev/report_gen_alg.jl +++ b/docs/src/dev/report_gen_alg.jl @@ -1,49 +1,9 @@ -#= -julia --project=docs -using Revise; include(joinpath("docs", "src", "dev", "report_gen_alg.jl")) - -Tested algs: +import ArgParse +import JLD2 - - `SSP22Heuns` - - `SSP33ShuOsher` - - `RK4` - - `ARK2GKC` - - `ARS111` - - `ARS121` - - `ARS122` - - `ARS222` - - `ARS232` - - `ARS233` - - `ARS343` - - `ARS443` - - `SSP222` - - `SSP322` - - `SSP332` - - `SSP333` - - `SSP433` - - `DBM453` - - `HOMMEM1` - - `IMKG232a` - - `IMKG232b` - - `IMKG242a` - - `IMKG242b` - - `IMKG243a` - - `IMKG252a` - - `IMKG252b` - - `IMKG253a` - - `IMKG253b` - - `IMKG254a` - - `IMKG254b` - - `IMKG254c` - - `IMKG342a` - - `IMKG343a` - - `SSPKnoth` -=# -using ClimaTimeSteppers -using Test -using InteractiveUtils: subtypes +include(joinpath(@__DIR__, "compute_convergence.jl")) +include(joinpath(pkgdir(ClimaTimeSteppers), "test", "problems.jl")) -import ArgParse function parse_commandline() s = ArgParse.ArgParseSettings() ArgParse.@add_arg_table s begin @@ -55,55 +15,57 @@ function parse_commandline() parsed_args = ArgParse.parse_args(ARGS, s) return parsed_args end +parsed_args = parse_commandline() +alg_str = parsed_args["alg"] +alg_name = getproperty(ClimaTimeSteppers, Symbol(alg_str))() -ENV["GKSwstype"] = "nul" # avoid displaying plots +convergence_results = Dict() -include(joinpath(@__DIR__, "compute_convergence.jl")) -include(joinpath(pkgdir(ClimaTimeSteppers), "test", "problems.jl")) -import ClimaTimeSteppers as CTS -all_subtypes(::Type{T}) where {T} = isabstracttype(T) ? vcat(all_subtypes.(subtypes(T))...) : [T] - -let # Convergence - parsed_args = parse_commandline() - alg_name = getproperty(CTS, Symbol(parsed_args["alg"]))() - # NOTE: Some imperfections in the convergence order for SSPKnoth are to be - # expected because we are not using the exact Jacobian +test_convergence!(convergence_results, alg_name, ark_analytic_nonlin_test_cts(Float64), 450) +test_convergence!(convergence_results, alg_name, ark_analytic_sys_test_cts(Float64), 200) +test_convergence!( + convergence_results, + alg_name, + ark_analytic_test_cts(Float64), + 25000; + super_convergence_algorithm_names = (ARS121(),), +) +test_convergence!(convergence_results, alg_name, onewaycouple_mri_test_cts(Float64), 9000; high_order_sample_shifts = 3) - export_convergence_results(alg_name, ark_analytic_nonlin_test_cts(Float64), 200) - export_convergence_results(alg_name, ark_analytic_sys_test_cts(Float64), 400) +test_convergence!( + convergence_results, + alg_name, + climacore_1Dheat_test_cts(Float64), + 200; + numerical_reference_algorithm_name = ARK548L2SA2(), +) +test_convergence!( + convergence_results, + alg_name, + climacore_1Dheat_test_implicit_cts(Float64), + 200; + high_order_sample_shifts = 2, + numerical_reference_algorithm_name = ARK548L2SA2(), + broken_tests = (SSP22Heuns(), SSP33ShuOsher(), RK4(), IMKG253b(), IMKG254a(), IMKG254b(), IMKG343a()), +) # This problem exceeds the CFL bounds of ERK methods, and it is unstable for several IMKG methods. +test_convergence!( + convergence_results, + alg_name, + climacore_2Dheat_test_cts(Float64), + 200; + numerical_reference_algorithm_name = ARK548L2SA2(), +) - export_convergence_results(alg_name, ark_analytic_test_cts(Float64), 40000; super_convergence = (ARS121(),)) - export_convergence_results(alg_name, onewaycouple_mri_test_cts(Float64), 10000; num_steps_scaling_factor = 5) - export_convergence_results( - alg_name, - climacore_1Dheat_test_cts(Float64), - 400; - num_steps_scaling_factor = 4, - numerical_reference_algorithm_name = ARS343(), - ) +mkpath("output") +JLD2.save_object(joinpath("output", "convergence_$alg_str.jld2"), convergence_results) - rosenbrock_schemes = filter(name -> name isa ClimaTimeSteppers.RosenbrockAlgorithmName, get_algorithm_names()) - if any(x -> alg_name isa typeof(x), rosenbrock_schemes) - export_convergence_results(alg_name, climacore_1Dheat_test_implicit_cts(Float64), 400) - end - export_convergence_results( - alg_name, - climacore_2Dheat_test_cts(Float64), - 600; - num_steps_scaling_factor = 4, - numerical_reference_algorithm_name = ARS343(), - ) - - # Unconstrained vs SSP results without limiters - parsed_args = parse_commandline() - alg_name = getproperty(CTS, Symbol(parsed_args["alg"]))() - ssp_schemes = all_subtypes(ClimaTimeSteppers.IMEXSSPRKAlgorithmName) - if any(x -> alg_name isa x, ssp_schemes) - test_unconstrained_vs_ssp_without_limiters(alg_name, ark_analytic_nonlin_test_cts(Float64), 200) - test_unconstrained_vs_ssp_without_limiters(alg_name, ark_analytic_sys_test_cts(Float64), 400) - test_unconstrained_vs_ssp_without_limiters(alg_name, ark_analytic_test_cts(Float64), 40000) - test_unconstrained_vs_ssp_without_limiters(alg_name, onewaycouple_mri_test_cts(Float64), 10000) - test_unconstrained_vs_ssp_without_limiters(alg_name, climacore_1Dheat_test_cts(Float64), 200) - test_unconstrained_vs_ssp_without_limiters(alg_name, climacore_2Dheat_test_cts(Float64), 200) - end +# Unconstrained vs SSP tests +if alg_name isa ClimaTimeSteppers.IMEXSSPRKAlgorithmName + test_unconstrained_vs_ssp_without_limiters(alg_name, ark_analytic_nonlin_test_cts(Float64), 450) + test_unconstrained_vs_ssp_without_limiters(alg_name, ark_analytic_sys_test_cts(Float64), 200) + test_unconstrained_vs_ssp_without_limiters(alg_name, ark_analytic_test_cts(Float64), 25000) + test_unconstrained_vs_ssp_without_limiters(alg_name, onewaycouple_mri_test_cts(Float64), 9000) + test_unconstrained_vs_ssp_without_limiters(alg_name, climacore_1Dheat_test_cts(Float64), 200) + test_unconstrained_vs_ssp_without_limiters(alg_name, climacore_1Dheat_test_implicit_cts(Float64), 200) + test_unconstrained_vs_ssp_without_limiters(alg_name, climacore_2Dheat_test_cts(Float64), 200) end diff --git a/docs/src/dev/summarize_convergence.jl b/docs/src/dev/summarize_convergence.jl index 63c7dd2d..6bd36805 100644 --- a/docs/src/dev/summarize_convergence.jl +++ b/docs/src/dev/summarize_convergence.jl @@ -1,151 +1,51 @@ -#= -julia --project=docs -using Revise; include(joinpath("docs", "src", "dev", "summarize_convergence.jl")) -=# - -using ClimaTimeSteppers +import Plots import JLD2 -using Test +using LaTeXStrings: latexstring using Printf: @sprintf -ENV["GKSwstype"] = "nul" # avoid displaying plots -using InteractiveUtils: subtypes +ENV["GKSwstype"] = "nul" # avoid displaying plots include(joinpath(@__DIR__, "compute_convergence.jl")) -include(joinpath(pkgdir(ClimaTimeSteppers), "test", "problems.jl")) - -all_subtypes(::Type{T}) where {T} = isabstracttype(T) ? vcat(all_subtypes.(subtypes(T))...) : [T] - -function _rprint_dict_structure(io::IO, x::T, pc, xname) where {T <: NamedTuple} - for pn in propertynames(x) - pc_full = (pc..., ".", pn) - xi = getproperty(x, pn) - _rprint_dict_structure(io, xi; pc = pc_full, xname) - end -end; - -function _rprint_dict_structure(io::IO, x::Dict; pc, xname) - for k in keys(x) - pc_full = (pc..., "[$k]") - xi = getindex(x, k) - _rprint_dict_structure(io, xi; pc = pc_full, xname) - end -end; - -function _rprint_dict_structure(io::IO, xi; pc, xname) - println(io, "$(xname * string(join(pc)))") -end - -_rprint_dict_structure(io::IO, x::T, xname) where {T <: Union{NamedTuple, Dict}} = - _rprint_dict_structure(io, x; pc = (), xname) -_rprint_dict_structure(x::T, xname) where {T <: Union{NamedTuple, Dict}} = _rprint_dict_structure(stdout, x, xname) - -""" - @rprint_dict_structure(::T) where {T <: Union{NamedTuple, Dict}} - -Recursively print keys of a (potentially nested) dict. -""" -macro rprint_dict_structure(x) - return :(_rprint_dict_structure(stdout, $(esc(x)), $(string(x)))) -end -function algorithm_names_by_availability(out_dict, test_name, algorithm_names_all, plot1_dts) - # @rprint_dict_structure out_dict # for debugging - algorithm_names = [] - for algorithm_name in algorithm_names_all - key2 = string(algorithm_name) - keep_alg = true - if !haskey(out_dict, key2) - keep_alg = false - @warn "out_dict has no key $key2" - else - if !haskey(out_dict[key2], "cur_avg_errs_dict") - keep_alg = false - @warn "out_dict[$key2] has no key cur_avg_errs_dict" - else - for plot1_dt in plot1_dts - if !haskey(out_dict[key2]["cur_avg_errs_dict"], plot1_dt) - keep_alg = false - @warn "out_dict[$key2][cur_avg_errs_dict] has no key plot1_dt" - end - end - end - end - keep_alg && push!(algorithm_names, algorithm_name) +float_str(x) = @sprintf "%.4f" x +pow_str(x) = "10^{$(@sprintf "%.1f" log10(x))}" +function si_str(x) + if isnan(x) || x in (0, Inf, -Inf) + return string(x) end - return algorithm_names + exponent = floor(Int, log10(x)) + mantissa = x / 10.0^exponent + return "$(float_str(mantissa)) \\times 10^{$exponent}" end -function summarize_convergence( - out_dict_test, - alg_name, - test_case, - num_steps; - num_steps_scaling_factor = 10, - order_confidence_percent = 99, - super_convergence = (), - numerical_reference_algorithm_name = nothing, - numerical_reference_num_steps = num_steps_scaling_factor^3 * num_steps, - full_history_algorithm_name = nothing, - average_function = array -> norm(array) / sqrt(length(array)), - average_function_str = "RMS", - only_endpoints = false, - verbose = false, -) - - title = "All Algorithms" - algorithm_names_all = get_algorithm_names() - - (; test_name, t_end, linear_implicit, analytic_sol) = test_case - - keep_alg = true - plot1_dts = t_end ./ round.(Int, num_steps .* num_steps_scaling_factor .^ (-1:0.5:1)) - algorithm_names = algorithm_names_by_availability(out_dict_test, test_name, algorithm_names_all, plot1_dts) - @show algorithm_names - - # out_dict_test = Dict() - # out_dict_test[key2] = Dict() - - prob = test_case.split_prob - FT = typeof(t_end) - default_dt = t_end / num_steps - - algorithm(algorithm_name::ClimaTimeSteppers.ERKAlgorithmName) = ExplicitAlgorithm(algorithm_name) - algorithm(algorithm_name::ClimaTimeSteppers.SSPKnoth) = - ClimaTimeSteppers.RosenbrockAlgorithm(ClimaTimeSteppers.tableau(ClimaTimeSteppers.SSPKnoth())) - algorithm(algorithm_name::ClimaTimeSteppers.IMEXARKAlgorithmName) = - IMEXAlgorithm(algorithm_name, NewtonsMethod(; max_iters = linear_implicit ? 1 : 2)) - - ref_sol = if isnothing(numerical_reference_algorithm_name) - analytic_sol - else - ref_alg = algorithm(numerical_reference_algorithm_name) - ref_alg_str = string(nameof(typeof(numerical_reference_algorithm_name))) - ref_dt = t_end / numerical_reference_num_steps - verbose && - @info "Generating numerical reference solution for $test_name with $ref_alg_str (dt = $ref_dt)..." - out_dict_test["numerical_ref_sol"] # solve(deepcopy(prob), ref_alg; dt = ref_dt, save_everystep = !only_endpoints) - end - - cur_avg_err(u, t, integrator) = average_function(abs.(u .- ref_sol(t))) - cur_avg_sol_and_err(u, t, integrator) = (average_function(u), average_function(abs.(u .- ref_sol(t)))) - - float_str(x) = @sprintf "%.4f" x - pow_str(x) = "10^{$(@sprintf "%.1f" log10(x))}" - function si_str(x) - if isnan(x) || x in (0, Inf, -Inf) - return string(x) - end - exponent = floor(Int, log10(x)) - mantissa = x / 10.0^exponent - return "$(float_str(mantissa)) \\times 10^{$exponent}" +function summarize_convergence(test_name, test_results, default_dt, numerical_reference_info) + rms_solution_str = "\\textrm{rms}\\_\\textrm{solution}" + rms_error_str = "\\textrm{rms}\\_\\textrm{error}" + average_rms_solution_str = "\\textrm{average}\\_$rms_solution_str" + average_rms_error_str = "\\textrm{average}\\_$rms_error_str" + default_dt_str = "\$dt = $(pow_str(default_dt))\$" + reference_str = + isnothing(numerical_reference_info) ? "Y_{analytic}(t)[\\textrm{index}]" : "Y_{ref}(t)[\\textrm{index}]" + rms_solution_def_str = "\\textrm{RMS}_{\\textrm{index}}\\{Y(t)[\\textrm{index}]\\}" + rms_error_def_str = "\\textrm{RMS}_{\\textrm{index}}\\{Y(t)[\\textrm{index}] - $reference_str\\}" + average_rms_solution_def_str = "\\textrm{RMS}_t\\{$rms_solution_str(t)\\}" + average_rms_error_def_str = "\\textrm{RMS}_t\\{$rms_error_str(t)\\}" + footnote = """ + Terminology: + \$$rms_solution_str(t) = $rms_solution_def_str\$ + \$$rms_error_str(t) = $rms_error_def_str\$ + \$$average_rms_solution_str = $average_rms_solution_def_str\$ + \$$average_rms_error_str = $average_rms_error_def_str\$""" + if !isnothing(numerical_reference_info) + (ref_alg_str, ref_dt, ref_average_rms_error) = numerical_reference_info + ref_error_def_str = "Y_{ref}(t)[\\textrm{index}] - Y_{analytic}(t)[\\textrm{index}]" + ref_average_rms_error_def_str = "\\textrm{RMS}_{t,\\,\\textrm{index}}\\{$ref_error_def_str\\}" + footnote = "$footnote\n\n\nThe numerical reference solution \$Y_{ref}\$ is \ + computed using $ref_alg_str with\n\$dt = $(pow_str(ref_dt)),\\ \ + \\textrm{and}\\ $ref_average_rms_error_def_str = \ + $(si_str(ref_average_rms_error))\$" end - net_avg_sol_str = "\\textrm{$average_function_str}\\_\\textrm{solution}" - net_avg_err_str = "\\textrm{$average_function_str}\\_\\textrm{error}" - cur_avg_sol_str = "\\textrm{current}\\_$net_avg_sol_str" - cur_avg_err_str = "\\textrm{current}\\_$net_avg_err_str" - linestyles = (:solid, :dash, :dot, :dashdot, :dashdotdot) marker_kwargs = (; markershape = :circle, markeralpha = 0.5, markerstrokewidth = 0) plot_kwargs = (; @@ -159,163 +59,110 @@ function summarize_convergence( bottommargin = 30Plots.px, ) + plot1_min = Inf + plot1_max = -Inf plot1 = Plots.plot(; title = "Convergence Orders", xaxis = (latexstring("dt"), :log10), - yaxis = (latexstring(net_avg_err_str), :log10), - legendtitle = "Convergence Order ($order_confidence_percent% CI)", + yaxis = (latexstring(average_rms_error_str), :log10), + legendtitle = "Convergence Order (99% CI)", plot_kwargs..., ) - plot2b_min = typemax(FT) - plot2b_max = typemin(FT) + plot2a_min = Inf + plot2a_max = -Inf plot2a = Plots.plot(; - title = latexstring("Solutions with \$dt = $(pow_str(default_dt))\$"), + title = latexstring("Solutions with $default_dt_str"), xaxis = (latexstring("t"),), - yaxis = (latexstring(cur_avg_sol_str),), - legendtitle = latexstring(net_avg_sol_str), + yaxis = (latexstring(rms_solution_str),), + legendtitle = latexstring(average_rms_solution_str), plot_kwargs..., ) + + plot2b_min = Inf + plot2b_max = -Inf plot2b = Plots.plot(; - title = latexstring("Errors with \$dt = $(pow_str(default_dt))\$"), + title = latexstring("Errors with $default_dt_str"), xaxis = (latexstring("t"),), - yaxis = (latexstring(cur_avg_err_str), :log10), - legendtitle = latexstring(net_avg_err_str), + yaxis = (latexstring(rms_error_str), :log10), + legendtitle = latexstring(average_rms_error_str), plot_kwargs..., ) - scb_cur_avg_err = make_saving_callback(cur_avg_err, prob.u0, t_end, nothing) - scb_cur_avg_sol_and_err = make_saving_callback(cur_avg_sol_and_err, prob.u0, t_end, nothing) + alg_strs = collect(keys(test_results)) + for alg_str in all_subtypes(ClimaTimeSteppers.AbstractAlgorithmName) + if !(string(alg_str) in alg_strs) + @warn "Convergence of $alg_str has not been tested with $test_name" + end + end + sort_by(alg_str) = ( + test_results[alg_str]["predicted_order"], + -round(test_results[alg_str]["average_rms_errors"][end]; sigdigits = 2), + alg_str, + ) # sort by predicted order, then by rounded error at large dt, then by name + sort!(alg_strs; by = sort_by) + for alg_str in alg_strs + predicted_order = test_results[alg_str]["predicted_order"] + predicted_super_convergence = test_results[alg_str]["predicted_super_convergence"] + sampled_dts = test_results[alg_str]["sampled_dts"] + average_rms_errors = test_results[alg_str]["average_rms_errors"] + default_dt_times = test_results[alg_str]["default_dt_times"] + default_dt_solutions = test_results[alg_str]["default_dt_solutions"] + default_dt_errors = test_results[alg_str]["default_dt_errors"] - for algorithm_name in algorithm_names - cur_avg_errs_dict = out_dict_test[string(algorithm_name)]["cur_avg_errs_dict"] - alg = algorithm(algorithm_name) - alg_str = string(nameof(typeof(algorithm_name))) - predicted_order = predicted_convergence_order(algorithm_name, prob.f) linestyle = linestyles[(predicted_order - 1) % length(linestyles) + 1] - - verbose && @info "Running $test_name with $alg_str..." - plot1_net_avg_errs = map(plot1_dts) do plot1_dt - cur_avg_errs = cur_avg_errs_dict[plot1_dt] - # solve( - # deepcopy(prob), - # alg; - # dt = plot1_dt, - # save_everystep = !only_endpoints, - # callback = scb_cur_avg_err, - # ).u - verbose && @info "RMS_error(dt = $plot1_dt) = $(average_function(cur_avg_errs))" - return average_function(cur_avg_errs) - end - order, order_uncertainty = convergence_order(plot1_dts, plot1_net_avg_errs, order_confidence_percent / 100) + order, order_uncertainty = convergence_order(sampled_dts, average_rms_errors, 0.99) order_str = "$(float_str(order)) \\pm $(float_str(order_uncertainty))" - if algorithm_name in super_convergence - predicted_order += 1 - plot1_label = "$alg_str: \$$order_str\\ \\ \\ \\textbf{\\textit{SC}}\$" - else - plot1_label = "$alg_str: \$$order_str\$" - end - verbose && @info "Order = $order ± $order_uncertainty" - if abs(order - predicted_order) > order_uncertainty - @warn "Predicted order outside error bars for $alg_str ($test_name)" - end - if order_uncertainty > predicted_order / 10 - @warn "Order uncertainty too large for $alg_str ($test_name)" - end - Plots.plot!(plot1, plot1_dts, plot1_net_avg_errs; label = latexstring(plot1_label), linestyle, marker_kwargs...) - # Remove all 0s from plot2_cur_avg_errs because they cannot be plotted on a - # logarithmic scale. Record the extrema of plot2_cur_avg_errs to set ylim. - plot2_data = out_dict_test[string(algorithm_name)]["plot2_data"] - # plot2_data = solve( - # deepcopy(prob), - # alg; - # dt = default_dt, - # save_everystep = !only_endpoints, - # callback = scb_cur_avg_sol_and_err, - # ) - plot2_ts = plot2_data.t - plot2_cur_avg_sols = first.(plot2_data.u) - plot2_cur_avg_errs = last.(plot2_data.u) - plot2b_min = min(plot2b_min, minimum(x -> x == 0 ? typemax(FT) : x, plot2_cur_avg_errs)) - plot2b_max = max(plot2b_max, maximum(plot2_cur_avg_errs)) - plot2_cur_avg_errs .= max.(plot2_cur_avg_errs, eps(FT(0))) - plot2_net_avg_sol = average_function(plot2_cur_avg_sols) - plot2_net_avg_err = average_function(plot2_cur_avg_errs) + # Ignore large values from algorithms without reasonable order estimates. + ignore_max = isnan(order_uncertainty) || order_uncertainty > 10 + + plot1_min = min(plot1_min, minimum(filter(!isnan, average_rms_errors); init = Inf)) + ignore_max || (plot1_max = max(plot1_max, maximum(average_rms_errors))) + plot1_label = + predicted_super_convergence ? "$alg_str: \$$order_str\\ \\ \\ \\textbf{\\textit{SC}}\$" : + "$alg_str: \$$order_str\$" + Plots.plot!( + plot1, + sampled_dts, + average_rms_errors; + label = latexstring(plot1_label), + linestyle, + marker_kwargs..., + ) + + plot2a_min = min(plot2a_min, minimum(filter(!isnan, default_dt_solutions))) + ignore_max || (plot2a_max = max(plot2a_max, maximum(default_dt_solutions))) Plots.plot!( plot2a, - plot2_ts, - plot2_cur_avg_sols; - label = latexstring("$alg_str: \$$(si_str(plot2_net_avg_sol))\$"), + default_dt_times, + default_dt_solutions; + label = latexstring("$alg_str: \$$(si_str(rms(default_dt_solutions)))\$"), linestyle, - (only_endpoints ? marker_kwargs : (;))..., ) + + # Remove all 0s from default_dt_errors because they cannot be plotted on + # a logarithmic scale. + plot2b_min = min(plot2b_min, minimum(filter(!iszero, filter(!isnan, default_dt_errors)))) + ignore_max || (plot2b_max = max(plot2b_max, maximum(default_dt_errors))) + default_dt_errors .= max.(default_dt_errors, eps(0.0)) Plots.plot!( plot2b, - plot2_ts, - plot2_cur_avg_errs; - label = latexstring("$alg_str: \$$(si_str(plot2_net_avg_err))\$"), + default_dt_times, + default_dt_errors; + label = latexstring("$alg_str: \$$(si_str(rms(default_dt_errors)))\$"), linestyle, - (only_endpoints ? marker_kwargs : (;))..., ) end + # Add a factor of 2 for padding in the log plots. + Plots.plot!(plot1; ylim = (plot1_min / 2, plot1_max * 2)) Plots.plot!(plot2b; ylim = (plot2b_min / 2, plot2b_max * 2)) - plots = (plot1, plot2a, plot2b) + # Add 3% of the range for padding in the linear plot. + plot2a_padding = (plot2a_max - plot2a_min) * 0.03 + Plots.plot!(plot2a; ylim = (plot2a_min - plot2a_padding, plot2a_max + plot2a_padding)) - if !isnothing(full_history_algorithm_name) - history_alg = algorithm(full_history_algorithm_name) - history_alg_name = string(nameof(typeof(full_history_algorithm_name))) - history_solve_results = out_dict_test[history_alg_name]["history_solve_results"] - # history_solve_results = solve( - # deepcopy(prob), - # history_alg; - # dt = default_dt, - # save_everystep = !only_endpoints, - # callback = make_saving_callback((u, t, integrator) -> u .- ref_sol(t), prob.u0, t_end, nothing), - # ) - history_array = hcat(history_solve_results.u...) - history_plot_title = "Errors for $history_alg_name with \$dt = $(pow_str(default_dt))\$" - history_plot = Plots.heatmap( - history_solve_results.t, - 1:size(history_array, 1), - history_array; - title = latexstring(history_plot_title), - xaxis = (latexstring("t"),), - yaxis = ("index",), - plot_kwargs..., - ) - plots = (plots..., history_plot) - end - - avg_def_str(val_str, var_str) = "\\textrm{$average_function_str}(\\{$val_str\\,\\forall\\,$var_str\\})" - ref_sol_def_str = - isnothing(numerical_reference_algorithm_name) ? "Y_{analytic}(t)[\\textrm{index}]" : - "Y_{ref}(t)[\\textrm{index}]" - cur_avg_sol_def_str = avg_def_str("Y(t)[\\textrm{index}]", "\\textrm{index}") - cur_avg_err_def_str = avg_def_str("|Y(t)[\\textrm{index}] - $ref_sol_def_str|", "\\textrm{index}") - net_avg_sol_def_str = avg_def_str("$cur_avg_sol_str(t)", "t") - net_avg_err_def_str = avg_def_str("$cur_avg_err_str(t)", "t") - footnote = """ - Terminology: - \$$cur_avg_sol_str(t) = $cur_avg_sol_def_str\$ - \$$cur_avg_err_str(t) = $cur_avg_err_def_str\$ - \$$net_avg_sol_str = $net_avg_sol_def_str\$ - \$$net_avg_err_str = $net_avg_err_def_str\$""" - if !isnothing(numerical_reference_algorithm_name) - ref_cur_avg_errs = map(ref_sol.u, ref_sol.t) do u, t - average_function(abs.(u .- analytic_sol(t))) - end - ref_net_avg_err = average_function(ref_cur_avg_errs) - ref_cur_avg_err_def_str = - avg_def_str("|Y_{ref}(t)[\\textrm{index}] - Y_{analytic}(t)[\\textrm{index}]|", "\\textrm{index}") - ref_net_avg_err_def_str = avg_def_str(ref_cur_avg_err_def_str, "t") - footnote = "$footnote\n\n\nNote: The \"reference solution\" \$Y_{ref}\$ was \ - computed using $ref_alg_str with\n\$dt = $(pow_str(ref_dt)),\\ \ - \\textrm{and}\\ $ref_net_avg_err_def_str = \ - $(si_str(ref_net_avg_err))\$" - end footnote_plot = Plots.plot(; title = latexstring(footnote), titlelocation = :left, @@ -325,52 +172,56 @@ function summarize_convergence( margin = 0Plots.px, bottommargin = 10Plots.px, ) - plots = (plots..., footnote_plot) - n_plots = length(plots) + all_plots = (plot1, plot2a, plot2b, footnote_plot) + n_plots = length(all_plots) heights₀ = repeat([1 / (n_plots - 1)], n_plots - 1) footnote_height = eps() heights = [heights₀ .- footnote_height / n_plots..., footnote_height] layout = Plots.grid(n_plots, 1; heights) - plot = Plots.plot( - plots...; - plot_title = "Analysis of $title for \"$test_name\"", + full_plot = Plots.plot( + all_plots...; + plot_title = "Convergence analysis for \"$test_name\"", layout = layout, fontfamily = "Computer Modern", ) mkpath("output") - file_suffix = lowercase(replace(test_name * ' ' * title, " " => "_")) - Plots.savefig(plot, joinpath("output", "convergence_$file_suffix.png")) + file_suffix = lowercase(replace(test_name, " " => "_")) + Plots.savefig(full_plot, joinpath("output", "convergence_$file_suffix.png")) end - -function merge_tests(v::Vector) - test_names = collect(Set(map(x -> first(keys(x)), v))) - return Dict(map(test_names) do tn - matched_dicts = filter(x -> first(keys(x)) == tn, v) - inner_dicts = map(x -> first(values(x)), matched_dicts) - inner_merged = merge(inner_dicts...) - Pair(tn, inner_merged) - end...) +convergence_results_filenames = filter(readdir("output"; join = true)) do filename + startswith(basename(filename), "convergence_") && endswith(filename, ".jld2") end - -let # Convergence - # NOTE: Some imperfections in the convergence order for SSPKnoth are to be - # expected because we are not using the exact Jacobian - vector_of_dicts = [] - files = readdir() - filter!(x -> endswith(x, ".jld2"), files) - filter!(x -> startswith(basename(x), "convergence_"), files) - for f in files - push!(vector_of_dicts, JLD2.load_object(f)) +all_convergence_results = map(JLD2.load_object, convergence_results_filenames) +representative_test_data = merge(all_convergence_results...) +test_names = collect(keys(representative_test_data)) +for test_name in test_names, convergence_results in all_convergence_results + test_name in keys(convergence_results) || continue + @assert convergence_results[test_name]["default_dt"] == representative_test_data[test_name]["default_dt"] + @assert convergence_results[test_name]["numerical_reference_info"] == + representative_test_data[test_name]["numerical_reference_info"] +end # Check that each test has a consistent default dt and reference solution +for test_name in test_names + all_alg_results = map(all_convergence_results) do convergence_results + test_name in keys(convergence_results) ? convergence_results[test_name]["all_alg_results"] : Dict() end + summarize_convergence( + test_name, + merge(all_alg_results...), + representative_test_data[test_name]["default_dt"], + representative_test_data[test_name]["numerical_reference_info"], + ) +end - out_dict = merge_tests(vector_of_dicts) - for test_name in keys(out_dict) - out_dict_test = out_dict[test_name] - args = out_dict_test["args"] - kwargs = out_dict_test["kwargs"] - summarize_convergence(out_dict_test, args...; kwargs...) - end +#= TO TEST CONVERGENCE OF ALL ALGORITHMS: +include("docs/src/dev/compute_convergence.jl") +for alg_name in all_subtypes(AbstractAlgorithmName) + empty!(ARGS) + push!(ARGS, "--alg", string(alg_name)) + @info "Testing convergence of $alg_name" + include("docs/src/dev/report_gen_alg.jl") end +include("docs/src/dev/summarize_convergence.jl") +=# diff --git a/docs/src/plotting_utils.jl b/docs/src/plotting_utils.jl deleted file mode 100644 index 00c0422c..00000000 --- a/docs/src/plotting_utils.jl +++ /dev/null @@ -1,190 +0,0 @@ -import Plots, Markdown -using ClimaCorePlots -using Distributions: quantile, TDist -using Printf: @sprintf -using LaTeXStrings: latexstring -using PrettyTables: pretty_table, ft_printf -import ClimaTimeSteppers as CTS -import DiffEqCallbacks - -""" - predicted_convergence_order(algorithm_name, ode_function) - -Return the predicted convergence order of the algorithm for the given ODE -function (assuming that the algorithm converges). -""" -function predicted_convergence_order(algorithm_name::AbstractAlgorithmName, ode_function::AbstractClimaODEFunction) - (imp_order, exp_order, combined_order) = imex_convergence_orders(algorithm_name) - has_imp = !isnothing(ode_function.T_imp!) - has_exp = CTS.has_T_exp(ode_function) - has_imp && !has_exp && return imp_order - !has_imp && has_exp && return exp_order - has_imp && has_exp && return combined_order - return 0 -end - -""" - imex_convergence_orders(algorithm_name) - -Return a tuple containing the expected convergence order of the tableau when -using only an implicit tendency, the order when using only an explicit tendency, -and the order when using both tendencies. -""" -function imex_convergence_orders end -imex_convergence_orders(::ARS111) = (1, 1, 1) -imex_convergence_orders(::ARS121) = (1, 1, 1) -imex_convergence_orders(::ARS122) = (2, 2, 2) -imex_convergence_orders(::ARS222) = (2, 2, 2) -imex_convergence_orders(::ARS232) = (2, 3, 2) -imex_convergence_orders(::ARS233) = (3, 3, 3) -imex_convergence_orders(::ARS343) = (3, 4, 3) -imex_convergence_orders(::ARS443) = (3, 3, 3) -imex_convergence_orders(::Union{IMKG232a, IMKG232b}) = (2, 2, 2) -imex_convergence_orders(::Union{IMKG242a, IMKG242b}) = (2, 4, 2) -imex_convergence_orders(::IMKG243a) = (2, 4, 2) -imex_convergence_orders(::Union{IMKG252a, IMKG252b}) = (2, 2, 2) -imex_convergence_orders(::Union{IMKG253a, IMKG253b}) = (2, 2, 2) -imex_convergence_orders(::Union{IMKG254a, IMKG254b, IMKG254c}) = (2, 2, 2) -imex_convergence_orders(::IMKG342a) = (3, 4, 3) -imex_convergence_orders(::IMKG343a) = (3, 4, 3) -imex_convergence_orders(::SSP222) = (2, 2, 2) -imex_convergence_orders(::SSP322) = (2, 2, 2) -imex_convergence_orders(::SSP332) = (2, 3, 2) -imex_convergence_orders(::SSP333) = (3, 3, 3) -imex_convergence_orders(::SSP433) = (3, 3, 3) -imex_convergence_orders(::DBM453) = (3, 3, 3) -imex_convergence_orders(::HOMMEM1) = (2, 3, 2) -imex_convergence_orders(::ARK2GKC) = (2, 2, 2) -imex_convergence_orders(::ARK437L2SA1) = (4, 4, 4) -imex_convergence_orders(::ARK548L2SA2) = (5, 5, 5) -imex_convergence_orders(::SSP22Heuns) = (2, 2, 2) -imex_convergence_orders(::SSP33ShuOsher) = (3, 3, 3) -imex_convergence_orders(::RK4) = (4, 4, 4) -# SSPKnoth is not really an IMEX method -imex_convergence_orders(::SSPKnoth) = (2, 2, 2) - -# Compute a confidence interval for the convergence order, returning the -# estimated convergence order and its uncertainty. -function convergence_order(dts, errs, confidence) - log_dts = log10.(dts) - log_errs = log10.(errs) - n_dts = length(dts) - - # slope and vertical intercept of linear regression (for each log_dt and - # log_err, log_err ≈ log_dt * order + log_err_of_dt1) - order, log_err_of_dt1 = hcat(log_dts, ones(n_dts)) \ log_errs - - # number of degrees of freedom of linear regression (number of data - # points minus number of fitted parameters) - n_dof = n_dts - 2 - - # critical value of Student's t-distribution for two-sided confidence - # interval - confidence_interval_t_value = -quantile(TDist(n_dof), (1 - confidence) / 2) - - # standard deviation of linear regression - regression_standard_deviation = sqrt(sum((log_errs .- (log_dts .* order .+ log_err_of_dt1)) .^ 2) / n_dof) - - # standard deviation of slope - order_standard_deviation = regression_standard_deviation / sqrt(sum((log_dts .- sum(log_dts) / n_dts) .^ 2)) - - # "uncertainty" in slope (half of width of confidence interval) - order_uncertainty = confidence_interval_t_value * order_standard_deviation - - return order, order_uncertainty -end - -function make_saving_callback(cb, u, t, integrator) - savevalType = typeof(cb(u, t, integrator)) - return DiffEqCallbacks.SavingCallback(cb, DiffEqCallbacks.SavedValues(typeof(t), savevalType)) -end - -# Generates the Table 1 from -# "Optimization-based limiters for the spectral element method" by Guba et al., -# and also plots the values used to generate the table. -function limiter_summary(sol_dicts, algorithm_names) - to_title(name) = titlecase(replace(string(name), '_' => ' ')) - table_rows = [] - mkpath("output") - for alg_name in algorithm_names - alg_str = string(nameof(typeof(alg_name))) - plots = [] - plot_kwargs = (; - clims = (0, 1), - color = :diverging_rainbow_bgymr_45_85_c67_n256, - colorbar = false, - guide = "", - margin = 10Plots.px, - ) - for use_limiter in (false, true), use_hyperdiffusion in (false, true) - solution = sol_dicts[alg_name, use_hyperdiffusion, use_limiter] - initial_q = solution[1].ρq ./ solution[1].ρ - final_q = solution[end].ρq ./ solution[end].ρ - names = propertynames(initial_q) - - if isempty(plots) - for name in names - push!(plots, Plots.plot(initial_q.:($name); plot_kwargs..., title = to_title(name))) - end - end - for name in names - push!(plots, Plots.plot(final_q.:($name); plot_kwargs..., title = "")) - end - - for name in names - ϕ₀ = initial_q.:($name) - ϕ = final_q.:($name) - Δϕ₀ = maximum(ϕ₀) - minimum(ϕ₀) - ϕ_error = ϕ .- ϕ₀ - table_row = [ - alg_str;; - string(use_limiter);; - string(use_hyperdiffusion);; - to_title(name);; - max(0, -(minimum(ϕ) - minimum(ϕ₀)) / Δϕ₀);; - max(0, (maximum(ϕ) - maximum(ϕ₀)) / Δϕ₀);; - map(p -> norm(ϕ_error, p) / norm(ϕ₀, p), (1, 2, Inf))... - ] - push!(table_rows, table_row) - end - end - colorbar_plot = Plots.scatter( - [0]; - plot_kwargs..., - colorbar = true, - framestyle = :none, - legend_position = :none, - margin = 0Plots.px, - markeralpha = 0, - zcolor = [0], - ) - plot = Plots.plot( - plots..., - colorbar_plot; - layout = (Plots.@layout [Plots.grid(5, 3) a{0.1w}]), - plot_title = "Tracer specific humidity for $alg_str (Initial, \ - Final, Final w/ Hyperdiffusion, Final w/ Limiter, \ - Final w/ Hyperdiffusion & Limiter)", - size = (1600, 2000), - ) - Plots.savefig(plot, joinpath("output", "limiter_summary_$alg_str.png")) - end - table = pretty_table( - vcat(table_rows...); - header = [ - "Algorithm", - "Limiter", - "Hyperdiffusion", - "Tracer Name", - "Max Undershoot", - "Max Overshoot", - "1-Norm Error", - "2-Norm Error", - "∞-Norm Error", - ], - crop = :none, - body_hlines = collect(3:3:(length(table_rows) - 1)), - formatters = ft_printf("%.4e"), - ) - println(table) -end diff --git a/ext/ClimaTimeSteppersBenchmarkToolsExt.jl b/ext/ClimaTimeSteppersBenchmarkToolsExt.jl index 0acebceb..e069ce51 100644 --- a/ext/ClimaTimeSteppersBenchmarkToolsExt.jl +++ b/ext/ClimaTimeSteppersBenchmarkToolsExt.jl @@ -108,7 +108,6 @@ function CTS.benchmark_step( @assert all(x -> x in allowed_names, only) "Allowed names in `only` are: $allowed_names" end - W = get_W(integrator) X = similar(u) Xlim = similar(u) @. X = u @@ -116,10 +115,15 @@ function CTS.benchmark_step( trials₀ = OrderedCollections.OrderedDict() kwargs = (; device, with_cu_prof, trace, crop, hcrop) #! format: off - maybe_push!(trials₀, "Wfact", wfact_fun(integrator), (W, u, p, dt, t), kwargs, only) - maybe_push!(trials₀, "ldiv!", LA.ldiv!, (X, W, u), kwargs, only) - maybe_push!(trials₀, "T_imp!", implicit_fun(integrator), implicit_args(integrator), kwargs, only) - maybe_push!(trials₀, "T_exp_T_lim!", remaining_fun(integrator), remaining_args(integrator), kwargs, only) + if !isnothing(implicit_fun(integrator)) + W = get_W(integrator) + maybe_push!(trials₀, "Wfact", wfact_fun(integrator), (W, u, p, dt, t), kwargs, only) + maybe_push!(trials₀, "ldiv!", LA.ldiv!, (X, W, u), kwargs, only) + maybe_push!(trials₀, "T_imp!", implicit_fun(integrator), implicit_args(integrator), kwargs, only) + end + if !isnothing(remaining_fun(integrator)) + maybe_push!(trials₀, "T_exp_T_lim!", remaining_fun(integrator), remaining_args(integrator), kwargs, only) + end maybe_push!(trials₀, "lim!", f.lim!, (Xlim, p, t, u), kwargs, only) maybe_push!(trials₀, "dss!", f.dss!, (u, p, t), kwargs, only) maybe_push!(trials₀, "cache!", f.cache!, (u, p, t), kwargs, only) diff --git a/test/problems.jl b/test/problems.jl index c2e034cf..ced6c0e8 100644 --- a/test/problems.jl +++ b/test/problems.jl @@ -439,8 +439,6 @@ function onewaycouple_mri_test_cts(::Type{FT}) where {FT} ) end -Wfact!(W, Y, p, dtγ, t) = nothing - """ climacore_2Dheat_test_cts(::Type{<:AbstractFloat}) @@ -491,21 +489,7 @@ function climacore_2Dheat_test_cts(::Type{FT}) where {FT} return state end - # we add implicit pieces here for inference analysis - T_lim! = (Yₜ, u, _, t) -> nothing - cache_imp! = (u, _, t) -> nothing - cache! = (u, _, t) -> nothing - - jacobian = ClimaCore.MatrixFields.FieldMatrix((@name(u), @name(u)) => FT(-1) * LinearAlgebra.I) - - T_imp! = SciMLBase.ODEFunction( - (Yₜ, u, _, t) -> (Yₜ .= 0); - jac_prototype = FieldMatrixWithSolver(jacobian, init_state), - Wfact = Wfact!, - tgrad = (∂Y∂t, Y, p, t) -> (∂Y∂t .= 0), - ) - - tendency_func = ClimaODEFunction(; T_exp!, T_imp!, dss!, cache_imp!, cache!) + tendency_func = ClimaODEFunction(; T_exp!, dss!) split_tendency_func = tendency_func make_prob(func) = ODEProblem(func, init_state, (FT(0), t_end), nothing) IntegratorTestCase( @@ -561,9 +545,9 @@ function climacore_1Dheat_test_cts(::Type{FT}) where {FT} end function climacore_1Dheat_test_implicit_cts(::Type{FT}) where {FT} - n_elem_z = 10000 + n_elem_z = 100 # this should be larger than the explicit problem's n_elem_z n_z = 1 - f_0 = FT(0.0) # denoted by f̂₀ above + f_0 = FT(0) # denoted by f̂₀ above Δλ = FT(1) # denoted by Δλ̂ above t_end = FT(0.1) # denoted by t̂ above @@ -577,23 +561,21 @@ function climacore_1Dheat_test_implicit_cts(::Type{FT}) where {FT} init_state = Fields.FieldVector(; u = φ_sin) - diverg = Operators.DivergenceC2F(; bottom = Operators.SetDivergence(FT(0)), top = Operators.SetDivergence(FT(0))) + div = Operators.DivergenceC2F(; bottom = Operators.SetDivergence(FT(0)), top = Operators.SetDivergence(FT(0))) grad = Operators.GradientF2C() + function T_imp_func!(tendency, state, _, t) + @. tendency.u = div(grad(state.u)) + f_0 * exp(-(λ + Δλ) * t) * φ_sin + end - diverg_matrix = ClimaCore.MatrixFields.operator_matrix(diverg) + div_matrix = ClimaCore.MatrixFields.operator_matrix(div) grad_matrix = ClimaCore.MatrixFields.operator_matrix(grad) - function Wfact(W, Y, p, dtγ, t) name = @name(u) # NOTE: We need MatrixFields.⋅, not LinearAlgebra.⋅ - @. W.matrix[name, name] = dtγ * diverg_matrix() ⋅ grad_matrix() - (LinearAlgebra.I,) + @. W.matrix[name, name] = dtγ * div_matrix() ⋅ grad_matrix() - (LinearAlgebra.I,) return nothing end - function T_imp_func!(tendency, state, _, t) - @. tendency.u = diverg.(grad.(state.u)) + f_0 * exp(-(λ + Δλ) * t) * φ_sin - end - function tgrad(∂Y∂t, state, _, t) @. ∂Y∂t.u = -f_0 * (λ + Δλ) * exp(-(λ + Δλ) * t) * φ_sin end