|
| 1 | +using Test |
| 2 | +import Printf |
| 3 | +import Plots |
| 4 | + |
| 5 | +""" |
| 6 | + test_algs( |
| 7 | + algs_name, |
| 8 | + algs, |
| 9 | + test_case, |
| 10 | + num_steps; |
| 11 | + save_every_n_steps = max(1, Int(fld(num_steps, 500))), |
| 12 | + no_increment_algs = (), |
| 13 | + ) |
| 14 | +
|
| 15 | +Check that all of the specified ODE algorithms have the predicted convergence |
| 16 | +order, and that the increment formulations give the same results as the tendency |
| 17 | +formulations. Generate plots that show the algorithms' solutions for the |
| 18 | +specified test case, the errors of these solutions, and the convergence of these |
| 19 | +errors with respect to `dt`. |
| 20 | +
|
| 21 | +# Arguments |
| 22 | +- `algs_name::String`: the name of the collection of algorithms |
| 23 | +- `tableaus`: an array of tableaus |
| 24 | +- `test_case::IntegratorTestCase`: the test case to use |
| 25 | +- `num_steps::Int`: the numerical solutions for the solution and error plots |
| 26 | + are computed with `dt = t_end / num_steps`, while the solutions for |
| 27 | + the convergence plots are computed with `dt`, `dt / sqrt(10)`, and |
| 28 | + `dt * sqrt(10)` |
| 29 | +- `save_every_n_steps::Int`: the solution and error plots show only show the |
| 30 | + values at every `n`-th step; the default value is such that 500-999 steps |
| 31 | + are plotted (unless there are fewer than 500 steps, in which case every |
| 32 | + step is plotted) |
| 33 | +""" |
| 34 | +function test_algs( |
| 35 | + algs_name, |
| 36 | + tableaus, |
| 37 | + test_case, |
| 38 | + num_steps; |
| 39 | + save_every_n_steps = Int(cld(num_steps, 500)), |
| 40 | + super_convergence = nothing, |
| 41 | +) |
| 42 | + (; test_name, linear_implicit, t_end, analytic_sol) = test_case |
| 43 | + FT = typeof(t_end) |
| 44 | + linestyles = (:solid, :dash, :dot, :dashdot, :dashdotdot) |
| 45 | + plot_kwargs = (; |
| 46 | + size = (1000, 600), |
| 47 | + margin = 4Plots.mm, |
| 48 | + titlelocation = :left, |
| 49 | + legend_position = :outerright, |
| 50 | + palette = :glasbey_bw_minc_20_maxl_70_n256, |
| 51 | + ) |
| 52 | + |
| 53 | + plot1_dt = t_end / num_steps |
| 54 | + plot1_saveat = 0:(plot1_dt * save_every_n_steps):t_end |
| 55 | + plot1a = Plots.plot(; |
| 56 | + title = "Solution Norms of $algs_name Methods for `$test_name` \ |
| 57 | + (with dt = 10^$(Printf.@sprintf "%.1f" log10(plot1_dt)))", |
| 58 | + xlabel = "t", |
| 59 | + ylabel = "Solution Norm: ||Y_computed||", |
| 60 | + plot_kwargs..., |
| 61 | + ) |
| 62 | + plot1b = Plots.plot(; |
| 63 | + title = "Solution Errors of $algs_name Methods for `$test_name` \ |
| 64 | + (with dt = 10^$(Printf.@sprintf "%.1f" log10(plot1_dt)))", |
| 65 | + xlabel = "t", |
| 66 | + ylabel = "Error Norm: ||Y_computed - Y_analytic||", |
| 67 | + yscale = :log10, |
| 68 | + plot_kwargs..., |
| 69 | + ) |
| 70 | + plot1b_ymin = typemax(FT) # dynamically set ylim because some errors are 0 |
| 71 | + plot1b_ymax = typemin(FT) |
| 72 | + |
| 73 | + t_end_string = t_end % 1 == 0 ? string(Int(t_end)) : Printf.@sprintf("%.2f", t_end) |
| 74 | + plot2_dts = [plot1_dt / sqrt(10), plot1_dt, plot1_dt * sqrt(10)] |
| 75 | + plot2 = Plots.plot(; |
| 76 | + title = "Convergence Orders of $algs_name Methods for `$test_name` \ |
| 77 | + (at t = $t_end_string)", |
| 78 | + xlabel = "dt", |
| 79 | + ylabel = "Error Norm: ||Y_computed - Y_analytic||", |
| 80 | + xscale = :log10, |
| 81 | + yscale = :log10, |
| 82 | + plot_kwargs..., |
| 83 | + ) |
| 84 | + |
| 85 | + analytic_sols = map(analytic_sol, plot1_saveat) |
| 86 | + analytic_end_sol = [analytic_sols[end]] |
| 87 | + |
| 88 | + for tab in tableaus |
| 89 | + (prob, alg) = problem_algo(test_case, tab) |
| 90 | + predicted_order = if super_convergence == tab |
| 91 | + CTS.theoretical_convergence_order(tab()) + 1 |
| 92 | + else |
| 93 | + CTS.theoretical_convergence_order(tab()) |
| 94 | + end |
| 95 | + linestyle = linestyles[(predicted_order - 1) % length(linestyles) + 1] |
| 96 | + alg_name = string(nameof(tab)) |
| 97 | + |
| 98 | + # Use tstops to fix saving issues due to machine precision (e.g. if the |
| 99 | + # integrator needs to save at t but it stops at t - eps(), it will skip |
| 100 | + # over saving at t, unless tstops forces it to round t - eps() to t). |
| 101 | + solve_args = (; dt = plot1_dt, saveat = plot1_saveat, tstops = plot1_saveat) |
| 102 | + tendency_sols = solve(deepcopy(prob), alg; solve_args...).u |
| 103 | + tendency_norms = @. norm(tendency_sols) |
| 104 | + tendency_errs = @. norm(tendency_sols - analytic_sols) |
| 105 | + min_err = minimum(x -> x == 0 ? typemax(FT) : x, tendency_errs) |
| 106 | + plot1b_ymin = min(plot1b_ymin, min_err) |
| 107 | + plot1b_ymax = max(plot1b_ymax, maximum(tendency_errs)) |
| 108 | + tendency_errs .= max.(tendency_errs, eps(FT(0))) # plotting 0 breaks the log scale |
| 109 | + Plots.plot!(plot1a, plot1_saveat, tendency_norms; label = alg_name, linestyle) |
| 110 | + Plots.plot!(plot1b, plot1_saveat, tendency_errs; label = alg_name, linestyle) |
| 111 | + |
| 112 | + tendency_end_sols = map(dt -> solve(deepcopy(prob), alg; dt).u[end], plot2_dts) |
| 113 | + tendency_end_errs = @. norm(tendency_end_sols - analytic_end_sol) |
| 114 | + _, computed_order = hcat(ones(length(plot2_dts)), log10.(plot2_dts)) \ log10.(tendency_end_errs) |
| 115 | + @test computed_order ≈ predicted_order rtol = 0.1 |
| 116 | + @info "(alg, computed, predicted) = $alg, $computed_order, $predicted_order" |
| 117 | + label = "$alg_name ($(Printf.@sprintf "%.3f" computed_order))" |
| 118 | + Plots.plot!(plot2, plot2_dts, tendency_end_errs; label, linestyle) |
| 119 | + end |
| 120 | + Plots.plot!(plot1b; ylim = (plot1b_ymin / 2, plot1b_ymax * 2)) |
| 121 | + |
| 122 | + mkpath("output") |
| 123 | + file_suffix = "$(test_name)_$(lowercase(replace(algs_name, " " => "_")))" |
| 124 | + Plots.savefig(plot1a, joinpath("output", "solutions_$(file_suffix).png")) |
| 125 | + Plots.savefig(plot1b, joinpath("output", "errors_$(file_suffix).png")) |
| 126 | + Plots.savefig(plot2, joinpath("output", "orders_$(file_suffix).png")) |
| 127 | +end |
![bors[bot]](https://avatars.githubusercontent.com/in/1847?v=4&size=40){kind=avatar}
0 commit comments