From c370f28a305ae6ba6167168e83b22ea3324e9bb1 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Mon, 25 Sep 2023 21:27:30 +0200 Subject: [PATCH] create some more figures some improvements related to plotting --- create_figures.jl | 83 ++++++++++++++--- docs/src/overview.md | 9 +- examples/bbm_bbm_1d/bbm_bbm_1d_basic.jl | 2 +- examples/bbm_bbm_1d/bbm_bbm_1d_relaxation.jl | 2 +- .../bbm_bbm_variable_bathymetry_1d_basic.jl | 2 +- ...bm_bbm_variable_bathymetry_1d_dingemans.jl | 2 +- ...m_bbm_variable_bathymetry_1d_relaxation.jl | 2 +- ...bm_variable_bathymetry_1d_well_balanced.jl | 2 +- .../svaerd_kalisch_1d_dingemans.jl | 2 +- .../svaerd_kalisch_1d_dingemans_relaxation.jl | 2 +- .../svaerd_kalisch_1d_dingemans_upwind.jl | 2 +- .../svaerd_kalisch_1d_well_balanced.jl | 2 +- plot_examples.jl | 38 ++++---- src/callbacks_step/analysis.jl | 18 ++-- src/equations/bbm_bbm_1d.jl | 2 +- .../bbm_bbm_variable_bathymetry_1d.jl | 2 +- src/equations/svaerd_kalisch_1d.jl | 2 +- src/visualization.jl | 91 ++++++++++--------- 18 files changed, 165 insertions(+), 100 deletions(-) diff --git a/create_figures.jl b/create_figures.jl index 272fa3c1..4f01368b 100644 --- a/create_figures.jl +++ b/create_figures.jl @@ -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) @@ -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", @@ -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", @@ -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() diff --git a/docs/src/overview.md b/docs/src/overview.md index 10b19023..459eb646 100644 --- a/docs/src/overview.md +++ b/docs/src/overview.md @@ -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) ``` @@ -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) ``` diff --git a/examples/bbm_bbm_1d/bbm_bbm_1d_basic.jl b/examples/bbm_bbm_1d/bbm_bbm_1d_basic.jl index 00ac9e29..f2640569 100644 --- a/examples/bbm_bbm_1d/bbm_bbm_1d_basic.jl +++ b/examples/bbm_bbm_1d/bbm_bbm_1d_basic.jl @@ -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 diff --git a/examples/bbm_bbm_1d/bbm_bbm_1d_relaxation.jl b/examples/bbm_bbm_1d/bbm_bbm_1d_relaxation.jl index b5c49100..842ea53d 100644 --- a/examples/bbm_bbm_1d/bbm_bbm_1d_relaxation.jl +++ b/examples/bbm_bbm_1d/bbm_bbm_1d_relaxation.jl @@ -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 diff --git a/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_basic.jl b/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_basic.jl index 84318d55..54530d6e 100644 --- a/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_basic.jl +++ b/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_basic.jl @@ -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) diff --git a/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_dingemans.jl b/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_dingemans.jl index 7500149c..20ed2c5d 100644 --- a/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_dingemans.jl +++ b/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_dingemans.jl @@ -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) diff --git a/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_relaxation.jl b/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_relaxation.jl index 5f2cd73f..6916bdd9 100644 --- a/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_relaxation.jl +++ b/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_relaxation.jl @@ -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) diff --git a/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_well_balanced.jl b/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_well_balanced.jl index c591185f..11fe5bf2 100644 --- a/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_well_balanced.jl +++ b/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_well_balanced.jl @@ -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) diff --git a/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans.jl b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans.jl index e1ffcd54..c09c5478 100644 --- a/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans.jl +++ b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans.jl @@ -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) diff --git a/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans_relaxation.jl b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans_relaxation.jl index 0bac8033..52d7e038 100644 --- a/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans_relaxation.jl +++ b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans_relaxation.jl @@ -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) diff --git a/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans_upwind.jl b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans_upwind.jl index f4d548cd..3f3366a1 100644 --- a/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans_upwind.jl +++ b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans_upwind.jl @@ -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 diff --git a/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_well_balanced.jl b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_well_balanced.jl index 80a8a963..2844c82a 100644 --- a/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_well_balanced.jl +++ b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_well_balanced.jl @@ -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) diff --git a/plot_examples.jl b/plot_examples.jl index 00c8e606..3fb4e6fe 100644 --- a/plot_examples.jl +++ b/plot_examples.jl @@ -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) @@ -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 @@ -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)) ############################################################################### @@ -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)) ############################################################################### @@ -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)) ############################################################################### @@ -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)) ############################################################################### @@ -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)) ############################################################################### @@ -124,7 +124,7 @@ 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)) ############################################################################### @@ -132,8 +132,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_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), @@ -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), @@ -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), @@ -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), @@ -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)) diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index 8ff8e6da..75c16533 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -329,12 +329,12 @@ function analyze(quantity::Union{typeof(energy_total_modified), typeof(entropy_m integrate_quantity(quantity, u_ode, semi) end -pretty_form_utf(::typeof(waterheight_total)) = "∑η" -pretty_form_utf(::typeof(velocity)) = "∑v" -pretty_form_utf(::typeof(momentum)) = "∑P" -pretty_form_utf(::typeof(discharge)) = "∑P" -pretty_form_utf(::typeof(entropy)) = "∑U" -pretty_form_utf(::typeof(energy_total)) = "∑e_total" -pretty_form_utf(::typeof(entropy_modified)) = "∑U_modified" -pretty_form_utf(::typeof(energy_total_modified)) = "∑e_modified" -pretty_form_utf(::typeof(lake_at_rest_error)) = "∑|η-η₀|" +pretty_form_utf(::typeof(waterheight_total)) = "∫η" +pretty_form_utf(::typeof(velocity)) = "∫v" +pretty_form_utf(::typeof(momentum)) = "∫P" +pretty_form_utf(::typeof(discharge)) = "∫P" +pretty_form_utf(::typeof(entropy)) = "∫U" +pretty_form_utf(::typeof(energy_total)) = "∫e_total" +pretty_form_utf(::typeof(entropy_modified)) = "∫U_modified" +pretty_form_utf(::typeof(energy_total_modified)) = "∫e_modified" +pretty_form_utf(::typeof(lake_at_rest_error)) = "∫|η-η₀|" diff --git a/src/equations/bbm_bbm_1d.jl b/src/equations/bbm_bbm_1d.jl index c0411f5a..30446360 100644 --- a/src/equations/bbm_bbm_1d.jl +++ b/src/equations/bbm_bbm_1d.jl @@ -27,7 +27,7 @@ function BBMBBMEquations1D(; gravity_constant, D = 1.0) BBMBBMEquations1D(gravity_constant, D) end -varnames(::typeof(prim2prim), ::BBMBBMEquations1D) = ("eta", "v") +varnames(::typeof(prim2prim), ::BBMBBMEquations1D) = ("η", "v") varnames(::typeof(prim2cons), ::BBMBBMEquations1D) = ("h", "hv") # TODO: Initial condition should not get a `mesh` diff --git a/src/equations/bbm_bbm_variable_bathymetry_1d.jl b/src/equations/bbm_bbm_variable_bathymetry_1d.jl index 44e0587b..784b9c9d 100644 --- a/src/equations/bbm_bbm_variable_bathymetry_1d.jl +++ b/src/equations/bbm_bbm_variable_bathymetry_1d.jl @@ -27,7 +27,7 @@ function BBMBBMVariableEquations1D(; gravity_constant, eta0 = 1.0) BBMBBMVariableEquations1D(gravity_constant, eta0) end -varnames(::typeof(prim2prim), ::BBMBBMVariableEquations1D) = ("eta", "v", "D") +varnames(::typeof(prim2prim), ::BBMBBMVariableEquations1D) = ("η", "v", "D") varnames(::typeof(prim2cons), ::BBMBBMVariableEquations1D) = ("h", "hv", "b") # TODO: Initial condition should not get a `mesh` diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl index 22437d9a..68e6f945 100644 --- a/src/equations/svaerd_kalisch_1d.jl +++ b/src/equations/svaerd_kalisch_1d.jl @@ -41,7 +41,7 @@ function SvaerdKalischEquations1D(; gravity_constant, eta0 = 1.0, alpha = 0.0, SvaerdKalischEquations1D(gravity_constant, eta0, alpha, beta, gamma) end -varnames(::typeof(prim2prim), ::SvaerdKalischEquations1D) = ("eta", "v", "D") +varnames(::typeof(prim2prim), ::SvaerdKalischEquations1D) = ("η", "v", "D") varnames(::typeof(prim2cons), ::SvaerdKalischEquations1D) = ("h", "hv", "b") # TODO: Initial condition should not get a `mesh` diff --git a/src/visualization.jl b/src/visualization.jl index 318be639..62567773 100644 --- a/src/visualization.jl +++ b/src/visualization.jl @@ -1,23 +1,19 @@ -struct PlotData{Ylim, Conversion} +struct PlotData{Conversion} semisol::Pair{<:Semidiscretization, <:ODESolution} plot_initial::Bool plot_bathymetry::Bool conversion::Conversion step::Integer - # call this yli because ylim, ylimits etc. are already occupied by Plots.jl, but we want a vector - yli::Ylim end -struct PlotDataOverTime{RealT, Ylim, Conversion} +struct PlotDataOverTime{RealT, Conversion} semisol::Pair{<:Semidiscretization, <:ODESolution} x::RealT conversion::Conversion - # call this yli because ylim, ylimits etc. are already occupied by Plots.jl, but we want a vector - yli::Ylim end @recipe function f(plotdata::PlotData) - @unpack semisol, plot_initial, plot_bathymetry, conversion, step, yli = plotdata + @unpack semisol, plot_initial, plot_bathymetry, conversion, step = plotdata semi, sol = semisol equations = semi.equations names = varnames(conversion, equations) @@ -29,13 +25,9 @@ end nsubplots -= 1 end - if isnothing(yli) - yli = fill(:auto, nsubplots) - end if step == -1 step = length(sol.t) end - @assert length(yli)==nsubplots "The vector yli must be as long as there are subplots" initial_condition = semi.initial_condition t = sol.t[step] @@ -60,7 +52,6 @@ end end plot_title --> "$(get_name(semi.equations)) at t = $(round(t, digits=5))" - size --> (1200, 800) layout := nsubplots for i in 1:nsubplots @@ -70,6 +61,7 @@ end if plot_initial == true @series begin subplot := i + linestyle := :solid label := "initial $(names[i])" grid(semi), q_exact[i, :] end @@ -78,10 +70,9 @@ end @series begin subplot := i label --> names[i] - xguide := "x" - yguide := names[i] - title := names[i] - ylim := yli[i] + xguide --> "x" + yguide --> names[i] + title --> names[i] grid(semi), data[i, :] end end @@ -95,7 +86,6 @@ end xguide := "x" yguide := names[1] title := names[1] - ylim := yli[1] color := :black grid(semi), bathy end @@ -103,7 +93,7 @@ end end @recipe function f(plotdata_over_time::PlotDataOverTime) - @unpack semisol, x, conversion, yli = plotdata_over_time + @unpack semisol, x, conversion = plotdata_over_time semi, sol = semisol equations = semi.equations names = varnames(conversion, equations) @@ -115,10 +105,6 @@ end nsubplots -= 1 end - if isnothing(yli) - yli = fill(:auto, nsubplots) - end - index = argmin(abs.(grid(semi) .- x)) solution = zeros(nvariables(semi), length(sol.t)) data = zeros(nvars, length(sol.t)) @@ -147,44 +133,66 @@ end xguide := "t" yguide := names[i] title := names[i] - ylim := yli[i] sol.t, data[i, :] end end end @recipe function f(semisol::Pair{<:Semidiscretization, <:ODESolution}; plot_initial = false, - plot_bathymetry = true, conversion = prim2prim, step = -1, yli = nothing) - PlotData(semisol, plot_initial, plot_bathymetry, conversion, step, yli) + plot_bathymetry = true, conversion = prim2prim, step = -1) + PlotData(semisol, plot_initial, plot_bathymetry, conversion, step) end @recipe function f(semi::Semidiscretization, sol::ODESolution; plot_initial = false, - plot_bathymetry = true, conversion = prim2prim, step = -1, yli = nothing) - PlotData(semi => sol, plot_initial, plot_bathymetry, conversion, step, yli) + plot_bathymetry = true, conversion = prim2prim, step = -1) + PlotData(semi => sol, plot_initial, plot_bathymetry, conversion, step) end @recipe function f(semisol::Pair{<:Semidiscretization, <:ODESolution}, x_value; - conversion = prim2prim, yli = nothing) - PlotDataOverTime(semisol, x_value, conversion, yli) + conversion = prim2prim) + PlotDataOverTime(semisol, x_value, conversion) end @recipe function f(semi::Semidiscretization, sol::ODESolution, x_value; - conversion = prim2prim, yli = nothing) - PlotDataOverTime(semi => sol, x_value, conversion, yli) + conversion = prim2prim) + PlotDataOverTime(semi => sol, x_value, conversion) +end + +function pretty(name) + if name == :waterheight_total + return "∫η" + elseif name == :velocity + return "∫v" + elseif name in [:discharge, :momentum] + return "∫P" + elseif name == :entropy + return "∫U" + elseif name == :l2_error + return "L² error" + elseif name == :linf_error + return "L∞ error" + elseif name == :conservation_error + return "∫|q_q₀|" + else + return string(name) + end end @recipe function f(cb::DiscreteCallback{Condition, Affect!}; what = (:integrals,), - label_extension = "") where {Condition, Affect! <: AnalysisCallback} + label_extension = "", + exclude = []) where {Condition, Affect! <: AnalysisCallback} t = tstops(cb) subplot = 1 - layout := length(what) + layout --> length(what) if :integrals in what ints = integrals(cb) - plot_title --> "change of invariants" - for (name, integral) in pairs(ints) + + for (i, (name, integral)) in enumerate(pairs(ints)) + name in exclude && continue @series begin - subplot := subplot - label := string(name) * " " * label_extension + subplot --> subplot + label := pretty(name) * " " * label_extension + title --> "change of invariants" xguide --> "t" yguide --> "change of invariants" t, integral .- integral[1] @@ -194,11 +202,12 @@ end end if :errors in what errs = errors(cb) - plot_title --> "errors" - for (name, err) in pairs(errs) + for (i, (name, err)) in enumerate(pairs(errs)) + name in exclude && continue @series begin - subplot := subplot - label := string(name) * " " * label_extension + subplot --> subplot + label := pretty(name) * " " * label_extension + title --> "errors" xguide --> "t" yguide --> "sum of errors" t, dropdims(sum(err, dims = 1), dims = 1)