Skip to content

Commit 277f3a9

Browse files
MarcoArtianovchuravy
authored andcommitted
add quick convergence analysis for TRBDF2 as DIRK
1 parent 0cd6d3f commit 277f3a9

File tree

1 file changed

+56
-0
lines changed

1 file changed

+56
-0
lines changed

examples/trixi_convergence.jl

Lines changed: 56 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,56 @@
1+
# # Using the an implicit solver based on Ariadne with Trixi.jl
2+
3+
using Trixi
4+
using Implicit
5+
using CairoMakie
6+
7+
8+
# Notes:
9+
# Must disable both Polyester and LoopVectorization for Enzyme to be able to differentiate Trixi.jl
10+
# Using https://github.com/trixi-framework/Trixi.jl/pull/2295
11+
#
12+
# LocalPreferences.jl
13+
# ```toml
14+
# [Trixi]
15+
# loop_vectorization = false
16+
# polyester = false
17+
# ```
18+
19+
@assert !Trixi._PREFERENCE_POLYESTER
20+
@assert !Trixi._PREFERENCE_LOOPVECTORIZATION
21+
22+
trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_basic.jl"), cfl = 1.25/2, sol = nothing, save_solution = nothing,
23+
initial_refinement_level = 6, advection_velocity = (1.0, 1.0));
24+
25+
ode = semidiscretize(semi, (0.0, 2.0))
26+
27+
###############################################################################
28+
# run the simulation
29+
30+
sol = solve(
31+
ode,
32+
# Implicit.RKImplicitEuler();
33+
# Implicit.KS22();
34+
# Implicit.C23();
35+
# Implicit.L33();
36+
Implicit.RKTRBDF2();
37+
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
38+
ode_default_options()..., callback = callbacks,
39+
# verbose=1,
40+
krylov_algo = :gmres,
41+
);
42+
43+
cfl_vec = [10.0, 5.0, 2.5, 1.25, 1.25/2]
44+
l2_vec = [1.07287321e-02, 2.69376612e-03, 6.75387935e-04, 1.68974528e-04, 4.22696680e-05]
45+
linf_vec = [1.51727157e-02, 3.80957386e-03, 9.55254494e-04, 2.38966657e-04, 5.97799432e-05]
46+
f = Figure()
47+
ax = Axis(f[1, 1], xlabel = "CFL", ylabel = "Error",
48+
xscale = log10, yscale = log10, title = "Convergence Plot")
49+
50+
lines!(ax, cfl_vec, l2_vec, label = "L2 norm")
51+
lines!(ax, cfl_vec, linf_vec, label = "L∞ norm")
52+
lines!(ax, cfl_vec, 1e-4.*cfl_vec.^2, label = "2nd order")
53+
54+
axislegend(ax; position = :rb)
55+
56+
f

0 commit comments

Comments
 (0)