Skip to content

Commit

Permalink
create some more figures
Browse files Browse the repository at this point in the history
some improvements related to plotting
  • Loading branch information
JoshuaLampert committed Sep 25, 2023
1 parent 55275a3 commit c370f28
Show file tree
Hide file tree
Showing 18 changed files with 165 additions and 100 deletions.
83 changes: 71 additions & 12 deletions create_figures.jl
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,7 @@ function fig_2()
plot(k, c_euler.(k) ./ c0, label = "Euler", ylim = ylim, xguide = L"k",
yguide = L"c/c_0", linewidth = linewidth, markershape = :circle,
markersize = markersize)
plot!(k_zoom, c_euler.(k_zoom) ./ c0, ylims = (0.54, 1.0),
plot!(k_zoom, c_euler.(k_zoom) ./ c0, ylim = (0.54, 1.0),
inset = bbox(0.35, 0.1, 0.35, 0.3), subplot = 2, legend = nothing,
linewidth = linewidth, markershape = :circle, markersize = markersize)

Expand Down Expand Up @@ -112,18 +112,72 @@ function fig_2()
savefig(joinpath(OUT, "dispersion_relations.pdf"))
end

const OUT_SOLITON = joinpath(OUT, "soliton")
ispath(OUT_SOLITON) || mkpath(OUT_SOLITON)

# Plot errors, change of invariants, and solution at final time for baseline and relaxation
function fig_3_4_5()
linewidth = 2

g = 9.81
D = 2.0
c = 5 / 2 * sqrt(g * D)
x_min = -35.0
x_max = 35.0
tspan = (0.0, 50 * (x_max - x_min) / c)
N = 512
accuracy_order = 8

# baseline
trixi_include(joinpath(examples_dir(), EXAMPLES_DIR_BBMBBM, "bbm_bbm_1d_basic.jl"),
gravity_constant = g, D = D, coordinates_min = x_min,
coordinates_max = x_max, tspan = tspan, N = N,
accuracy_order = accuracy_order)
p1 = plot(analysis_callback, title = "", label_extension = "baseline", style = :auto,
linewidth = linewidth, layout = 2, subplot = 1)
p2 = plot(analysis_callback, title = "", what = (:errors,),
label_extension = "baseline", linestyle = :dash, linewidth = linewidth,
ylabel = L"\Vert\eta - \eta_{ana}\Vert_2 + \Vert v - v_{ana}\Vert_2",
exclude = [:conservation_error])
p3 = plot(semi => sol, label = "baseline", plot_initial = true, linestyle = :dash,
linewidth = linewidth, plot_title = "", title = "")

# relaxation
trixi_include(joinpath(examples_dir(), EXAMPLES_DIR_BBMBBM, "bbm_bbm_1d_relaxation.jl"),
gravity_constant = g, D = D, coordinates_min = x_min,
coordinates_max = x_max, tspan = tspan, N = N,
accuracy_order = accuracy_order)
plot!(p1, analysis_callback, title = "", label_extension = "relaxation", style = :auto,
linewidth = linewidth, subplot = 2)
plot!(p2, analysis_callback, title = "", what = (:errors,),
label_extension = "relaxation", linestyle = :dot, linewidth = linewidth,
ylabel = L"\Vert\eta - \eta_{ana}\Vert_2 + \Vert v - v_{ana}\Vert_2",
exclude = [:conservation_error])
plot!(p3, semi => sol, plot_bathymetry = false, label = "relaxation", linestyle = :dot,
linewidth = linewidth, plot_title = "", title = "", color = :green)

savefig(p1, joinpath(OUT_SOLITON, "invariants.pdf"))
savefig(p2, joinpath(OUT_SOLITON, "errors.pdf"))
savefig(p3, joinpath(OUT_SOLITON, "solution.pdf"))
end

# Plot convergence orders for baseline and relaxation
function fig_3()
function fig_6()
tspan = (0.0, 10.0)
accuracy_orders = [2, 4, 6, 8]
styles = [:dash, :dot, :dashdot, :dashdotdot]
iters = [4, 4, 4, 3]
initial_Ns = [128, 128, 128, 128]

all_Ns = minimum(initial_Ns) * 2 .^ (0:(maximum(iters) - 1))

plot([], label = :none, xscale = :log2, yscale = :log10, xticks = all_Ns, xlabel = "N",
ylabel = L"||\eta - \eta_{ana}||_2 + ||v - v_{ana}||_2")
linewidth = 2
markersize = 5
markershapes = [:circle, :star5, :star8, :rtriangle]
plot(label = :none, xscale = :log2, yscale = :log10, xlabel = "N", ylim = (1e-5, 1e2),
ylabel = L"\Vert\eta - \eta_{ana}\Vert_2 + \Vert v - v_{ana}\Vert_2",
legend = :bottomleft, layout = 2)

# left subplot: baseline
for i in 1:length(accuracy_orders)
Ns = initial_Ns[i] * 2 .^ (0:(iters[i] - 1))
_, errormatrix = convergence_test("examples/bbm_bbm_1d/bbm_bbm_1d_basic.jl",
Expand All @@ -133,12 +187,13 @@ function fig_3()
l2_err = sum(errormatrix[:l2], dims = 2)
eocs = log.(l2_err[2:end] ./ l2_err[1:(end - 1)]) ./ log(0.5)
eoc_mean = round(sum(eocs) / length(eocs), digits = 2)
plot!(Ns, l2_err, style = styles[i], label = "p = $accuracy_order, EOC: $eoc_mean")
plot!(Ns, l2_err, label = "p = $accuracy_order, EOC: $eoc_mean",
markershape = markershapes[i], linewidth = linewidth, markersize = markersize,
subplot = 1)
end
savefig(joinpath(OUT, "orders.pdf"))
xticks!(all_Ns, string.(all_Ns), subplot = 1)

plot([], label = :none, xscale = :log2, yscale = :log10, xticks = all_Ns, xlabel = "N",
ylabel = L"||\eta - \eta_{ana}||_2 + ||v - v_{ana}||_2")
# right subplot: relaxation
for i in 1:length(accuracy_orders)
Ns = initial_Ns[i] * 2 .^ (0:(iters[i] - 1))
_, errormatrix = convergence_test("examples/bbm_bbm_1d/bbm_bbm_1d_relaxation.jl",
Expand All @@ -148,11 +203,15 @@ function fig_3()
l2_err = sum(errormatrix[:l2], dims = 2)
eocs = log.(l2_err[2:end] ./ l2_err[1:(end - 1)]) ./ log(0.5)
eoc_mean = round(sum(eocs) / length(eocs), digits = 2)
plot!(Ns, l2_err, style = styles[i], label = "p = $accuracy_order, EOC: $eoc_mean")
plot!(Ns, l2_err, label = "p = $accuracy_order, EOC: $eoc_mean",
markershape = markershapes[i], linewidth = linewidth, markersize = markersize,
subplot = 2)
end
savefig(joinpath(OUT, "orders_relaxation.pdf"))
xticks!(all_Ns, string.(all_Ns), subplot = 2)
savefig(joinpath(OUT_SOLITON, "orders.pdf"))
end

fig_1()
fig_2()
fig_3()
fig_3_4_5()
fig_6()
9 changes: 3 additions & 6 deletions docs/src/overview.md
Original file line number Diff line number Diff line change
Expand Up @@ -105,16 +105,13 @@ savefig("shoaling_solution.png")

![](shoaling_solution.png)

By default, this will plot the bathymetry, but not the initial (analytical) solution. You can adjust this by passing the boolean values `plot_bathymetry` (if `true` always plot to first subplot) and `plot_initial`. You can also provide a `conversion` function that converts the solution. A conversion function should take the values of the primitive variables `q` at one node, and the `equations` as input and should return an `SVector` of any length as output. For a user defined conversion function, there should also exist a function `varnames(conversion, equations)` that returns a `Tuple` of the variable names used for labelling. The conversion function can, e.g., be [`prim2cons`](@ref) or [`waterheight_total`](@ref) if one only wants to plot the total waterheight. The resulting plot will have one subplot for each of the returned variables of the conversion variable. By default, the conversion function is just [`prim2prim`](@ref), i.e. the identity. The limits of the y-axis can be set by the keyword argument `yli`, which should be given as a vector of `Tuple`s with the same length as there are subplots.

!!! note "Setting y limits"
`ylim` and similar arguments are already occupied by Plots.jl and do not work with different y limits for multiple subplots.
By default, this will plot the bathymetry, but not the initial (analytical) solution. You can adjust this by passing the boolean values `plot_bathymetry` (if `true` always plot to first subplot) and `plot_initial`. You can also provide a `conversion` function that converts the solution. A conversion function should take the values of the primitive variables `q` at one node, and the `equations` as input and should return an `SVector` of any length as output. For a user defined conversion function, there should also exist a function `varnames(conversion, equations)` that returns a `Tuple` of the variable names used for labelling. The conversion function can, e.g., be [`prim2cons`](@ref) or [`waterheight_total`](@ref) if one only wants to plot the total waterheight. The resulting plot will have one subplot for each of the returned variables of the conversion variable. By default, the conversion function is just [`prim2prim`](@ref), i.e. the identity.

Plotting an animation over time can, e.g., be done by the following command, which uses `step` to plot the solution at a specific time step.

```@example overview
anim = @animate for step in 1:length(sol.u)
plot(semi => sol, plot_initial = true, conversion = waterheight_total, step = step, xlim = (-50, 20), yli = [(-0.8, 0.1)])
plot(semi => sol, plot_initial = true, conversion = waterheight_total, step = step, xlim = (-50, 20), ylims = [(-0.8, 0.1)])
end
gif(anim, "shoaling_solution.gif", fps = 25)
```
Expand Down Expand Up @@ -213,7 +210,7 @@ callbacks = CallbackSet(relaxation_callback, analysis_callback)
sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7,
save_everystep = false, callback = callbacks, saveat = saveat)
anim = @animate for step in 1:length(sol.u)
plot(semi => sol, plot_initial = true, conversion = waterheight_total, step = step, xlim = (-50, 20), yli = [(-0.8, 0.1)])
plot(semi => sol, plot_initial = true, conversion = waterheight_total, step = step, xlim = (-50, 20), ylims = [(-0.8, 0.1)])
end
gif(anim, "shoaling_solution_dg.gif", fps = 25)
```
Expand Down
2 changes: 1 addition & 1 deletion examples/bbm_bbm_1d/bbm_bbm_1d_basic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ boundary_conditions = boundary_condition_periodic
coordinates_min = -35.0
coordinates_max = 35.0
N = 512
mesh = Mesh1D(coordinates_min, coordinates_max, N + 1)
mesh = Mesh1D(coordinates_min, coordinates_max, N)

# create solver with periodic SBP operators of accuracy order 4
accuracy_order = 4
Expand Down
2 changes: 1 addition & 1 deletion examples/bbm_bbm_1d/bbm_bbm_1d_relaxation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ boundary_conditions = boundary_condition_periodic
coordinates_min = -35.0
coordinates_max = 35.0
N = 512
mesh = Mesh1D(coordinates_min, coordinates_max, N + 1)
mesh = Mesh1D(coordinates_min, coordinates_max, N)

# create solver with periodic SBP operators of accuracy order 4
accuracy_order = 4
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ boundary_conditions = boundary_condition_periodic
coordinates_min = -35.0
coordinates_max = 35.0
N = 512
mesh = Mesh1D(coordinates_min, coordinates_max, N + 1)
mesh = Mesh1D(coordinates_min, coordinates_max, N)

# create solver with periodic SBP operators of accuracy order 4
solver = Solver(mesh, 4)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ boundary_conditions = boundary_condition_periodic
coordinates_min = -138.0
coordinates_max = 46.0
N = 512
mesh = Mesh1D(coordinates_min, coordinates_max, N + 1)
mesh = Mesh1D(coordinates_min, coordinates_max, N)

# create solver with periodic SBP operators of accuracy order 4
solver = Solver(mesh, 4)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ boundary_conditions = boundary_condition_periodic
coordinates_min = -1.0
coordinates_max = 1.0
N = 512
mesh = Mesh1D(coordinates_min, coordinates_max, N + 1)
mesh = Mesh1D(coordinates_min, coordinates_max, N)

# create solver with periodic SBP operators
D1 = periodic_derivative_operator(1, 4, mesh.xmin, mesh.xmax, mesh.N)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ boundary_conditions = boundary_condition_periodic
coordinates_min = -1.0
coordinates_max = 1.0
N = 512
mesh = Mesh1D(coordinates_min, coordinates_max, N + 1)
mesh = Mesh1D(coordinates_min, coordinates_max, N)

# create solver with periodic SBP operators
D1 = periodic_derivative_operator(1, 4, mesh.xmin, mesh.xmax, mesh.N)
Expand Down
2 changes: 1 addition & 1 deletion examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans.jl
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ boundary_conditions = boundary_condition_periodic
coordinates_min = -138.0
coordinates_max = 46.0
N = 512
mesh = Mesh1D(coordinates_min, coordinates_max, N + 1)
mesh = Mesh1D(coordinates_min, coordinates_max, N)

# create solver with periodic SBP operators of accuracy order 4
solver = Solver(mesh, 4)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ boundary_conditions = boundary_condition_periodic
coordinates_min = -138.0
coordinates_max = 46.0
N = 512
mesh = Mesh1D(coordinates_min, coordinates_max, N + 1)
mesh = Mesh1D(coordinates_min, coordinates_max, N)

# create solver with periodic SBP operators of accuracy order 4
solver = Solver(mesh, 4)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ boundary_conditions = boundary_condition_periodic
coordinates_min = -138.0
coordinates_max = 46.0
N = 512
mesh = Mesh1D(coordinates_min, coordinates_max, N + 1)
mesh = Mesh1D(coordinates_min, coordinates_max, N)

# create solver with periodic SBP operators of accuracy order 4
accuracy_order = 4
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@ boundary_conditions = boundary_condition_periodic
coordinates_min = -1.0
coordinates_max = 1.0
N = 200
mesh = Mesh1D(coordinates_min, coordinates_max, N + 1)
mesh = Mesh1D(coordinates_min, coordinates_max, N)

# create solver with periodic SBP operators of accuracy order 4
solver = Solver(mesh, 4)
Expand Down
38 changes: 19 additions & 19 deletions plot_examples.jl
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ macro plot_example(filename, args...)

# Plot solution
anim = @animate for step in 1:length(sol.u)
plot(semi => sol, plot_initial = true, step = step, yli = $ylims_gif)
plot(semi => sol, plot_initial = true, step = step, ylims = $ylims_gif)
end
gif(anim, joinpath(outdir, "solution.gif"), fps = 25)

Expand All @@ -39,7 +39,7 @@ macro plot_example(filename, args...)
# Plot at different x values over time
@assert size($x_values) == size($tlims)
for (i, x) in enumerate($x_values)
plot(semi => sol, x, xlim = $tlims[i], yli = $ylims_x)
plot(semi => sol, x, xlim = $tlims[i], ylims = $ylims_x)
savefig(joinpath(outdir, "solution_at_x_" * string(x) * ".pdf"))
end
end
Expand All @@ -66,19 +66,19 @@ const EXAMPLES_DIR_SVAERD_KALISCH = "svaerd_kalisch_1d"
# Travelling wave solution for one-dimensional BBM-BBM equations with periodic boundary conditions
# using periodic SBP operators
@plot_example(joinpath(EXAMPLES_DIR_BBMBBM, "bbm_bbm_1d_basic.jl"),
ylims_gif=[(-8, 4), :auto], tspan=(0.0, 50.0))
ylims_gif=[(-8, 4) :auto], tspan=(0.0, 50.0))

###############################################################################
# Travelling wave solution for one-dimensional BBM-BBM equations with periodic boundary conditions
# using discontinuously coupled Legendre SBP operators
@plot_example(joinpath(EXAMPLES_DIR_BBMBBM, "bbm_bbm_1d_dg.jl"),
ylims_gif=[(-4, 2), :auto])
ylims_gif=[(-4, 2) :auto])

###############################################################################
# Travelling wave solution for one-dimensional BBM-BBM equations with periodic boundary conditions
# using periodic SBP operators and relaxation, is energy-conservative
@plot_example(joinpath(EXAMPLES_DIR_BBMBBM, "bbm_bbm_1d_relaxation.jl"),
ylims_gif=[(-8, 4), (-10, 30)],
ylims_gif=[(-8, 4) (-10, 30)],
tspan=(0.0, 30.0))

###############################################################################
Expand All @@ -87,7 +87,7 @@ const EXAMPLES_DIR_SVAERD_KALISCH = "svaerd_kalisch_1d"
# as a constant. Should give the same result as "bbm_bbm_1d_basic.jl"
@plot_example(joinpath(EXAMPLES_DIR_BBMBBM_VARIABLE,
"bbm_bbm_variable_bathymetry_1d_basic.jl"),
ylims_gif=[(-8, 4), :auto],
ylims_gif=[(-8, 4) :auto],
tspan=(0.0, 50.0))

###############################################################################
Expand All @@ -96,7 +96,7 @@ const EXAMPLES_DIR_SVAERD_KALISCH = "svaerd_kalisch_1d"
# is energy-conservative. Uses periodic finite difference SBP operators.
@plot_example(joinpath(EXAMPLES_DIR_BBMBBM_VARIABLE,
"bbm_bbm_variable_bathymetry_1d_relaxation.jl"),
ylims_gif=[(-1.5, 6.0), (-10.0, 10.0)],
ylims_gif=[(-1.5, 6.0) (-10.0, 10.0)],
tspan=(0.0, 10.0))

###############################################################################
Expand All @@ -105,7 +105,7 @@ const EXAMPLES_DIR_SVAERD_KALISCH = "svaerd_kalisch_1d"
# is energy-conservative. Uses upwind discontinuously coupled SBP operators.
@plot_elixir(joinpath(EXAMPLES_DIR_BBMBBM_VARIABLE,
"bbm_bbm_variable_bathymetry_1d_dg_upwind_relaxation.jl"),
ylims_gif=[(-1.5, 6.0), (-10.0, 10.0)],
ylims_gif=[(-1.5, 6.0) (-10.0, 10.0)],
tspan=(0.0, 10.0))

###############################################################################
Expand All @@ -114,7 +114,7 @@ const EXAMPLES_DIR_SVAERD_KALISCH = "svaerd_kalisch_1d"
# is energy-conservative. Uses periodic finite difference discontinuously coupled SBP operators.
@plot_example(joinpath(EXAMPLES_DIR_BBMBBM_VARIABLE,
"bbm_bbm_variable_bathymetry_1d_upwind_relaxation.jl"),
ylims_gif=[(-1.5, 6.0), (-10.0, 10.0)],
ylims_gif=[(-1.5, 6.0) (-10.0, 10.0)],
tspan=(0.0, 10.0))

###############################################################################
Expand All @@ -124,16 +124,16 @@ const EXAMPLES_DIR_SVAERD_KALISCH = "svaerd_kalisch_1d"
# (exactly) constant in time.
@plot_example(joinpath(EXAMPLES_DIR_BBMBBM_VARIABLE,
"bbm_bbm_variable_bathymetry_1d_well_balanced.jl"),
ylims_gif=[(2.0 - 1e-3, 2.0 + 1e-3), (-1e-3, 1e-3)],
ylims_gif=[(2.0 - 1e-3, 2.0 + 1e-3) (-1e-3, 1e-3)],
tspan=(0.0, 10.0))

###############################################################################
# One-dimensional BBM-BBM equations with initial condition that models
# a wave make. This setup comes from experiments by W. M. Dingemans.
@plot_example(joinpath(EXAMPLES_DIR_BBMBBM_VARIABLE,
"bbm_bbm_variable_bathymetry_1d_dingemans.jl"),
ylims_gif=[(-0.1, 0.9), (-0.3, 0.3)],
ylims_x=[:auto, :auto],
ylims_gif=[(-0.1, 0.9) (-0.3, 0.3)],
ylims_x=[:auto :auto],
x_values=[3.04, 9.44, 20.04, 26.04, 30.44, 37.04],
tlims=[
(15.0, 45.0),
Expand All @@ -150,8 +150,8 @@ const EXAMPLES_DIR_SVAERD_KALISCH = "svaerd_kalisch_1d"
# a wave make. This setup comes from experiments by W. M. Dingemans.
@plot_example(joinpath(EXAMPLES_DIR_SVAERD_KALISCH,
"svaerd_kalisch_1d_dingemans.jl"),
ylims_gif=[(-0.1, 0.9), (-0.3, 0.3)],
ylims_x=[:auto, :auto],
ylims_gif=[(-0.1, 0.9) (-0.3, 0.3)],
ylims_x=[:auto :auto],
x_values=[3.04, 9.44, 20.04, 26.04, 30.44, 37.04],
tlims=[
(15.0, 45.0),
Expand All @@ -168,8 +168,8 @@ const EXAMPLES_DIR_SVAERD_KALISCH = "svaerd_kalisch_1d"
# a wave make. This setup comes from experiments by W. M. Dingemans.
@plot_example(joinpath(EXAMPLES_DIR_SVAERD_KALISCH,
"svaerd_kalisch_1d_dingemans_upwind.jl"),
ylims_gif=[(-0.1, 0.9), (-0.3, 0.3)],
ylims_x=[:auto, :auto],
ylims_gif=[(-0.1, 0.9) (-0.3, 0.3)],
ylims_x=[:auto :auto],
x_values=[3.04, 9.44, 20.04, 26.04, 30.44, 37.04],
tlims=[
(15.0, 45.0),
Expand All @@ -187,8 +187,8 @@ const EXAMPLES_DIR_SVAERD_KALISCH = "svaerd_kalisch_1d"
# to preserve the modified entropy.
@plot_example(joinpath(EXAMPLES_DIR_SVAERD_KALISCH,
"svaerd_kalisch_1d_dingemans_relaxation.jl"),
ylims_gif=[(-0.1, 0.9), (-0.3, 0.3)],
ylims_x=[:auto, :auto],
ylims_gif=[(-0.1, 0.9) (-0.3, 0.3)],
ylims_x=[:auto :auto],
x_values=[3.04, 9.44, 20.04, 26.04, 30.44, 37.04],
tlims=[
(15.0, 45.0),
Expand All @@ -207,5 +207,5 @@ const EXAMPLES_DIR_SVAERD_KALISCH = "svaerd_kalisch_1d"
# (exactly) constant in time.
@plot_example(joinpath(EXAMPLES_DIR_SVAERD_KALISCH,
"svaerd_kalisch_1d_well_balanced.jl"),
ylims=[(2.0 - 1e-3, 2.0 + 1e-3), (-1e-3, 1e-3)],
ylims=[(2.0 - 1e-3, 2.0 + 1e-3) (-1e-3, 1e-3)],
tspan=(0.0, 10.0))
Loading

0 comments on commit c370f28

Please sign in to comment.