From daeb6a721e05f49d4cb89ac572b6e05b71c383b7 Mon Sep 17 00:00:00 2001 From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Date: Sat, 30 Sep 2023 14:01:59 +0200 Subject: [PATCH] More figures (#50) * add some more figures * format * create more figures * several minor changes * more figures * fix typo --- .github/workflows/Format-check.yml | 2 +- create_figures.jl | 256 ------ docs/src/overview.md | 12 +- docs/src/reproduce.md | 8 +- examples/bbm_bbm_1d/bbm_bbm_1d_dg.jl | 10 +- .../bbm_bbm_variable_bathymetry_1d_basic.jl | 3 +- ...able_bathymetry_1d_dg_upwind_relaxation.jl | 20 +- ...bm_bbm_variable_bathymetry_1d_dingemans.jl | 5 +- ...m_bbm_variable_bathymetry_1d_relaxation.jl | 11 +- ...ariable_bathymetry_1d_upwind_relaxation.jl | 8 +- ...bm_variable_bathymetry_1d_well_balanced.jl | 7 +- .../svaerd_kalisch_1d_dingemans.jl | 5 +- .../svaerd_kalisch_1d_dingemans_relaxation.jl | 3 +- .../svaerd_kalisch_1d_dingemans_upwind.jl | 2 +- .../svaerd_kalisch_1d_well_balanced.jl | 5 +- src/DispersiveShallowWater.jl | 4 +- src/equations/bbm_bbm_1d.jl | 11 +- .../bbm_bbm_variable_bathymetry_1d.jl | 33 +- src/equations/equations.jl | 2 +- src/equations/svaerd_kalisch_1d.jl | 24 +- src/util.jl | 2 +- src/visualization.jl | 33 +- test/test_bbm_bbm_variable_bathymetry_1d.jl | 18 +- test/test_svaerd_kalisch_1d.jl | 30 +- test/test_unit.jl | 10 +- visualization/bbm_bbm_1d_widestencil.jl | 43 + .../bbm_bbm_variable_bathymetry_1d_upwind.jl | 45 ++ ..._bbm_variable_bathymetry_1d_widestencil.jl | 43 + visualization/create_figures.jl | 733 ++++++++++++++++++ .../elixir_shallowwater_1d_dingemans.jl | 59 ++ .../plot_examples.jl | 2 +- 31 files changed, 1065 insertions(+), 384 deletions(-) delete mode 100644 create_figures.jl create mode 100644 visualization/bbm_bbm_1d_widestencil.jl create mode 100644 visualization/bbm_bbm_variable_bathymetry_1d_upwind.jl create mode 100644 visualization/bbm_bbm_variable_bathymetry_1d_widestencil.jl create mode 100644 visualization/create_figures.jl create mode 100644 visualization/elixir_shallowwater_1d_dingemans.jl rename plot_examples.jl => visualization/plot_examples.jl (99%) diff --git a/.github/workflows/Format-check.yml b/.github/workflows/Format-check.yml index 9edf4ae6..aa6e808e 100644 --- a/.github/workflows/Format-check.yml +++ b/.github/workflows/Format-check.yml @@ -25,7 +25,7 @@ jobs: - name: Install JuliaFormatter and format run: | julia -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter"))' - julia -e 'using JuliaFormatter; format(["src", "test", "examples", "create_figures.jl", "plot_examples.jl"], verbose = true)' + julia -e 'using JuliaFormatter; format(["src", "test", "examples", "visualization"], verbose = true)' - name: Format check run: | julia -e ' diff --git a/create_figures.jl b/create_figures.jl deleted file mode 100644 index dd6c497a..00000000 --- a/create_figures.jl +++ /dev/null @@ -1,256 +0,0 @@ -using DispersiveShallowWater -using Plots -using LaTeXStrings - -const OUT = "out/" -ispath(OUT) || mkpath(OUT) -const EXAMPLES_DIR_BBMBBM = "bbm_bbm_1d" - -# Plot of bathymetry and waterheight -function fig_1() - L = 1.0 - n = 100 - x = LinRange(0.0, L, n) - fontsize = 20 - - # just pick some function for b and eta that look nice - H = 1.012 - - b(x) = x * cos.(3 * pi * x) + H - plot(x, b, color = :gray, fill = (0, 0.8, :gray), fillstyle = :/, linewidth = 3, - legend = nothing, ticks = nothing, border = :none) - - eta(x) = x / (x^2 + 1) * sin(2 * pi * x) + H + 1.5 - plot!(x, eta, color = :blue, fill = (b.(x), 0.4, :blue), linewidth = 3) - - x1 = 0.2 - plot!([x1, x1], [b(x1), eta(x1)], line = (Plots.Arrow(:open, :both, 2.5, 2.0), :black), - annotation = (x1 - 0.08, (eta(x1) + b(x1)) / 2, text(L"h(t, x)", fontsize)), - linewidth = 2) - x2 = 0.4 - plot!([x2, x2], [0.0, b(x2)], line = (Plots.Arrow(:open, :both, 2.5, 2.0), :black), - annotation = (x2 + 0.06, b(x2) / 2, text(L"b(x)", fontsize)), linewidth = 2) - x3 = 0.8 - plot!([x3, x3], [0.0, eta(x3)], line = (Plots.Arrow(:open, :both, 2.5, 2.0), :black), - annotation = (x3 - 0.08, eta(x3) / 2, text(L"\eta(t, x)", fontsize)), - linewidth = 2) - - savefig(joinpath(OUT, "bathymetry.pdf")) -end - -# Plot of diserpersion relations -function fig_2() - linewidth = 2 - markersize = 5 - - h0 = 1.0 - g = 1.0 - c0 = sqrt(g * h0) - - k = 0.01:0.5:(8 * pi) - k_zoom = 0.01:0.3:pi - ylim = (0.0, 1.1) - - omega_euler(k) = sqrt(g * k) * sqrt(tanh(h0 * k)) - c_euler(k) = omega_euler(k) / k - 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, 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, - framestyle = :box) - - function plot_dispersion_relation(omega, label, markershape) - c(k) = omega(k) / k - plot!(k, c.(k) ./ c0, label = label, linewidth = linewidth, - markershape = markershape, markersize = markersize) - plot!(k_zoom, c.(k_zoom) ./ c0, subplot = 2, linewidth = linewidth, - markershape = markershape, markersize = markersize) - end - - omega_bbmbbm_(k, d0) = sqrt(g * h0) * k / (1 + 1 / 6 * (d0 * k)^2) - omega_bbmbbm(k) = omega_bbmbbm_(k, h0) - plot_dispersion_relation(omega_bbmbbm, "BBM-BBM", :cross) - - alpha_set1 = -1 / 3 * c0 * h0^2 - beta_set1 = 0.0 * h0^3 - gamma_set1 = 0.0 * c0 * h0^3 - - alpha_set2 = 0.0004040404040404049 * c0 * h0^2 - beta_set2 = 0.49292929292929294 * h0^3 - gamma_set2 = 0.15707070707070708 * c0 * h0^3 - - alpha_set3 = 0.0 * c0 * h0^2 - beta_set3 = 0.27946992481203003 * h0^3 - gamma_set3 = 0.0521077694235589 * c0 * h0^3 - - alpha_set4 = 0.0 * c0 * h0^2 - beta_set4 = 0.2308939393939394 * h0^3 - gamma_set4 = 0.04034343434343434 * c0 * h0^3 - - function char_equation(alpha, beta, gamma, k) - a = (1 + beta / h0 * k^2) - b = (-alpha - beta * alpha / h0 * k^2 - gamma / h0) * k^3 - c = -g * h0 * k^2 + gamma * alpha / h0 * k^6 - omega1 = (-b + sqrt(b^2 - 4 * a * c)) / (2 * a) - # omega2 = (-b - sqrt(b^2 - 4*a*c))/(2*a) - return omega1 - end - - omega_set1(k) = char_equation(alpha_set1, beta_set1, gamma_set1, k) - plot_dispersion_relation(omega_set1, "S.-K. set 1", :rtriangle) - - omega_set2(k) = char_equation(alpha_set2, beta_set2, gamma_set2, k) - plot_dispersion_relation(omega_set2, "S.-K. set 2", :star5) - - omega_set3(k) = char_equation(alpha_set3, beta_set3, gamma_set3, k) - plot_dispersion_relation(omega_set3, "S.-K. set 3", :star8) - - omega_set4(k) = char_equation(alpha_set4, beta_set4, gamma_set4, k) - plot_dispersion_relation(omega_set4, "S.-K. set 4", :diamond) - - # Plot box - plot!([0.0, pi], [0.54, 0.54], color = :black, label = :none) - plot!([0.0, pi], [1.0, 1.0], color = :black, label = :none) - plot!([0.0, 0.0], [0.54, 1.0], color = :black, label = :none) - plot!([pi, pi], [0.54, 1.0], color = :black, label = :none) - - # Plot connecting lines - plot!([pi, 6.8], [0.54, 0.629], color = :black, label = :none) - plot!([pi, 6.8], [1, 1.01], color = :black, label = :none) - - savefig(joinpath(OUT, "dispersion_relations.pdf")) -end - -const OUT_SOLITON = joinpath(OUT, "soliton") -ispath(OUT_SOLITON) || mkpath(OUT_SOLITON) - -# Plot convergence orders for baseline and relaxation -function fig_3() - tspan = (0.0, 10.0) - accuracy_orders = [2, 4, 6, 8] - iters = [4, 4, 4, 3] - initial_Ns = [128, 128, 128, 128] - - all_Ns = minimum(initial_Ns) * 2 .^ (0:(maximum(iters) - 1)) - - 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", - iters[i]; N = initial_Ns[i], tspan = tspan, - accuracy_order = accuracy_orders[i]) - # Use sum over all L^2-errors for all variables, i.e. ||η - η_ana||_2 + ||v - v_ana||_2 - 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, label = "p = $accuracy_order, EOC: $eoc_mean", - markershape = markershapes[i], linewidth = linewidth, markersize = markersize, - subplot = 1) - end - xticks!(all_Ns, string.(all_Ns), subplot = 1) - - # 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", - iters[i]; N = initial_Ns[i], tspan = tspan, - accuracy_order = accuracy_orders[i]) - # Use sum over all L^2-errors for all variables, i.e. ||η - η_ana||_2 + ||v - v_ana||_2 - 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, label = "p = $accuracy_order, EOC: $eoc_mean", - markershape = markershapes[i], linewidth = linewidth, markersize = markersize, - subplot = 2) - end - xticks!(all_Ns, string.(all_Ns), subplot = 2) - savefig(joinpath(OUT_SOLITON, "orders.pdf")) -end - -# Plot errors, change of invariants, and solution at final time for baseline and relaxation -function fig_4_5_6() - 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 = "", - ylims = [(-8, 3) (-1, 40)]) - x = DispersiveShallowWater.grid(semi) - q = DispersiveShallowWater.wrap_array(sol.u[end], semi) - plot!(p3, x, view(q, 1, :), inset = (1, bbox(0.11, 0.6, 0.35, 0.32)), subplot = 3, - xlim = (-20, -10), - ylim = (-0.05, 0.05), legend = nothing, linewidth = linewidth, linestyle = :dash, - color = 2, - tickfontsize = 5, yticks = [0.04, 0.0, -0.04], xticks = [-20, -15, -10], - framestyle = :box) - q_exact = DispersiveShallowWater.wrap_array(DispersiveShallowWater.compute_coefficients(initial_condition, - tspan[2], - semi), - semi) - plot!(p3, x, view(q_exact, 1, :), subplot = 3, legend = nothing, linewidth = linewidth, - linestyle = :dot, color = 1) - # Plot box - plot!(p3, [-20, -10], [-0.1, -0.1], color = :black, label = :none) - plot!(p3, [-20, -10], [0.1, 0.1], color = :black, label = :none) - plot!(p3, [-20, -20], [-0.1, 0.1], color = :black, label = :none) - plot!(p3, [-10, -10], [-0.1, 0.1], color = :black, label = :none) - - # Plot connecting lines - plot!(p3, [-20, -29], [-0.1, -3.6], color = :black, label = :none) - plot!(p3, [-10, -3.15], [-0.1, -3.6], color = :black, label = :none) - - # 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 = 3) - x = DispersiveShallowWater.grid(semi) - q = DispersiveShallowWater.wrap_array(sol.u[end], semi) - plot!(p3, x, view(q, 1, :), subplot = 3, legend = nothing, linewidth = linewidth, - linestyle = :dot, color = 3) - - savefig(p1, joinpath(OUT_SOLITON, "invariants.pdf")) - savefig(p2, joinpath(OUT_SOLITON, "errors.pdf")) - savefig(p3, joinpath(OUT_SOLITON, "solution.pdf")) -end - -fig_1() -fig_2() -fig_3() -fig_4_5_6() diff --git a/docs/src/overview.md b/docs/src/overview.md index f9e1b030..efcef62a 100644 --- a/docs/src/overview.md +++ b/docs/src/overview.md @@ -173,16 +173,16 @@ As an example, let us create a semidiscretization based on discontinuous Galerki ```@example overview mesh = Mesh1D(coordinates_min, coordinates_max, N) accuracy_order = 4 -Dop = legendre_derivative_operator(-1.0, 1.0, accuracy_order) -sbp_mesh = UniformPeriodicMesh1D(mesh.xmin, mesh.xmax, div(mesh.N, accuracy_order)) +D_legendre = legendre_derivative_operator(-1.0, 1.0, accuracy_order) +uniform_mesh = UniformPeriodicMesh1D(mesh.xmin, mesh.xmax, div(mesh.N, accuracy_order)) ``` Upwind DG operators in negative, central and positive operators can be obtained by `couple_discontinuously` ```@example overview -central = couple_discontinuously(Dop, sbp_mesh) -minus = couple_discontinuously(Dop, sbp_mesh, Val(:minus)) -plus = couple_discontinuously(Dop, sbp_mesh, Val(:plus)) +central = couple_discontinuously(D_legendre, uniform_mesh) +minus = couple_discontinuously(D_legendre, uniform_mesh, Val(:minus)) +plus = couple_discontinuously(D_legendre, uniform_mesh, Val(:plus)) D1 = PeriodicUpwindOperators(minus, central, plus) ``` @@ -224,7 +224,7 @@ For more details see also the [documentation of SummationByPartsOperators.jl](ht Some more examples sorted by the simulated equations can be found in the [examples/](https://github.com/JoshuaLampert/DispersiveShallowWater.jl/tree/main/examples) subdirectory. Especially, in [examples/svaerd\_kalisch\_1d/](https://github.com/JoshuaLampert/DispersiveShallowWater.jl/tree/main/examples/svaerd_kalisch_1d) you can find Julia scripts that solve the [`SvaerdKalischEquations1D`](@ref) that were not covered in this tutorial. The same steps as described above, however, apply in the same way to these equations. Attention must be paid for these equations because they do not conserve the classical total entropy ``\mathcal E``, but a modified entropy ``\hat{\mathcal E}``, available as [`entropy_modified`](@ref). -More examples, especially focussing on plotting, can be found in the scripts [create_figures.jl](https://github.com/JoshuaLampert/DispersiveShallowWater.jl/blob/main/create_figures.jl) and [plot_examples.jl](https://github.com/JoshuaLampert/DispersiveShallowWater.jl/blob/main/plot_examples.jl). +More examples, especially focussing on plotting, can be found in the scripts [create_figures.jl](https://github.com/JoshuaLampert/DispersiveShallowWater.jl/blob/main/visualization/create_figures.jl) and [plot_examples.jl](https://github.com/JoshuaLampert/DispersiveShallowWater.jl/blob/main/visualization/plot_examples.jl). ## References diff --git a/docs/src/reproduce.md b/docs/src/reproduce.md index ac409dcb..235be748 100644 --- a/docs/src/reproduce.md +++ b/docs/src/reproduce.md @@ -6,7 +6,13 @@ In order to reproduce all figures used in the master thesis "Structure-Preservin julia> using DispersiveShallowWater julia> include(DispersiveShallowWater.path_create_figures()) ``` -Executing this script may take a while. It will generate a folder `out/` with certain subfolders containing the figures. If you want to modify the plots or only produce a subset of plots, you can download the script [`create_figures.jl`](https://github.com/JoshuaLampert/DispersiveShallowWater.jl/blob/main/create_figures.jl), modify it accordingly and run it by: +Note that for one figure [Trixi.jl](https://github.com/trixi-framework/Trixi.jl) is needed, so download Trixi.jl first: + +```julia +julia> using Pkg +julia> Pkg.add("Trixi") +``` +Executing the script may take a while. It will generate a folder `out/` with certain subfolders containing the figures. If you want to modify the plots or only produce a subset of plots, you can download the script [`create_figures.jl`](https://github.com/JoshuaLampert/DispersiveShallowWater.jl/blob/main/visualization/create_figures.jl), modify it accordingly and run it by: ```julia julia> include("create_figures.jl") diff --git a/examples/bbm_bbm_1d/bbm_bbm_1d_dg.jl b/examples/bbm_bbm_1d/bbm_bbm_1d_dg.jl index d0fa282a..7d7d6543 100644 --- a/examples/bbm_bbm_1d/bbm_bbm_1d_dg.jl +++ b/examples/bbm_bbm_1d/bbm_bbm_1d_dg.jl @@ -22,11 +22,11 @@ mesh = Mesh1D(coordinates_min, coordinates_max, N) # create solver p = 3 # N needs to be divisible by p + 1 -Dop = legendre_derivative_operator(-1.0, 1.0, p + 1) -sbp_mesh = UniformPeriodicMesh1D(coordinates_min, coordinates_max, div(N, p + 1)) -D1 = couple_discontinuously(Dop, sbp_mesh) -D_pl = couple_discontinuously(Dop, sbp_mesh, Val(:plus)) -D_min = couple_discontinuously(Dop, sbp_mesh, Val(:minus)) +D_legendre = legendre_derivative_operator(-1.0, 1.0, p + 1) +uniform_mesh = UniformPeriodicMesh1D(coordinates_min, coordinates_max, div(N, p + 1)) +D1 = couple_discontinuously(D_legendre, uniform_mesh) +D_pl = couple_discontinuously(D_legendre, uniform_mesh, Val(:plus)) +D_min = couple_discontinuously(D_legendre, uniform_mesh, Val(:minus)) D2 = sparse(D_pl) * sparse(D_min) solver = Solver(D1, D2) 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 54530d6e..1e02e94e 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 @@ -17,7 +17,8 @@ N = 512 mesh = Mesh1D(coordinates_min, coordinates_max, N) # create solver with periodic SBP operators of accuracy order 4 -solver = Solver(mesh, 4) +accuracy_order = 4 +solver = Solver(mesh, accuracy_order) # semidiscretization holds all the necessary data structures for the spatial discretization semi = Semidiscretization(mesh, equations, initial_condition, solver, diff --git a/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_dg_upwind_relaxation.jl b/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_dg_upwind_relaxation.jl index 5409edbf..b2d19d38 100644 --- a/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_dg_upwind_relaxation.jl +++ b/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_dg_upwind_relaxation.jl @@ -9,23 +9,23 @@ using SparseArrays: sparse equations = BBMBBMVariableEquations1D(gravity_constant = 9.81) -# initial_condition_variable_bathymetry needs periodic boundary conditions -initial_condition = initial_condition_sin_bathymetry +# initial_condition_convergence_test needs periodic boundary conditions +initial_condition = initial_condition_convergence_test boundary_conditions = boundary_condition_periodic # create homogeneous mesh -coordinates_min = -1.0 -coordinates_max = 1.0 +coordinates_min = -35.0 +coordinates_max = 35.0 N = 512 mesh = Mesh1D(coordinates_min, coordinates_max, N) # create solver -accuracy_order = 4 -Dop = legendre_derivative_operator(-1.0, 1.0, accuracy_order) -sbp_mesh = UniformPeriodicMesh1D(mesh.xmin, mesh.xmax, div(mesh.N, accuracy_order)) -central = couple_discontinuously(Dop, sbp_mesh) -minus = couple_discontinuously(Dop, sbp_mesh, Val(:minus)) -plus = couple_discontinuously(Dop, sbp_mesh, Val(:plus)) +p = 3 # N needs to be divisible by p + 1 +D_legendre = legendre_derivative_operator(-1.0, 1.0, p + 1) +uniform_mesh = UniformPeriodicMesh1D(mesh.xmin, mesh.xmax, div(mesh.N, p + 1)) +central = couple_discontinuously(D_legendre, uniform_mesh) +minus = couple_discontinuously(D_legendre, uniform_mesh, Val(:minus)) +plus = couple_discontinuously(D_legendre, uniform_mesh, Val(:plus)) D1 = PeriodicUpwindOperators(minus, central, plus) D2 = sparse(plus) * sparse(minus) solver = Solver(D1, D2) 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 20ed2c5d..8a0a25a8 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 @@ -4,7 +4,7 @@ using DispersiveShallowWater ############################################################################### # Semidiscretization of the BBM-BBM equations -equations = BBMBBMVariableEquations1D(gravity_constant = 9.81) +equations = BBMBBMVariableEquations1D(gravity_constant = 9.81, eta0 = 0.8) initial_condition = initial_condition_dingemans boundary_conditions = boundary_condition_periodic @@ -16,7 +16,8 @@ N = 512 mesh = Mesh1D(coordinates_min, coordinates_max, N) # create solver with periodic SBP operators of accuracy order 4 -solver = Solver(mesh, 4) +accuracy_order = 4 +solver = Solver(mesh, accuracy_order) # semidiscretization holds all the necessary data structures for the spatial discretization semi = Semidiscretization(mesh, equations, initial_condition, solver, 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 6916bdd9..25fd604c 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 @@ -8,18 +8,19 @@ using SparseArrays: sparse equations = BBMBBMVariableEquations1D(gravity_constant = 9.81) -# initial_condition_variable_bathymetry needs periodic boundary conditions -initial_condition = initial_condition_sin_bathymetry +# initial_condition_convergence_test needs periodic boundary conditions +initial_condition = initial_condition_convergence_test boundary_conditions = boundary_condition_periodic # create homogeneous mesh -coordinates_min = -1.0 -coordinates_max = 1.0 +coordinates_min = -35.0 +coordinates_max = 35.0 N = 512 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) +accuracy_order = 4 +D1 = periodic_derivative_operator(1, accuracy_order, mesh.xmin, mesh.xmax, mesh.N) D2 = sparse(D1)^2 solver = Solver(D1, D2) diff --git a/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_upwind_relaxation.jl b/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_upwind_relaxation.jl index 26458a49..b645b5ce 100644 --- a/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_upwind_relaxation.jl +++ b/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_upwind_relaxation.jl @@ -8,13 +8,13 @@ using SparseArrays: sparse equations = BBMBBMVariableEquations1D(gravity_constant = 9.81) -# initial_condition_variable_bathymetry needs periodic boundary conditions -initial_condition = initial_condition_sin_bathymetry +# initial_condition_convergence_test needs periodic boundary conditions +initial_condition = initial_condition_convergence_test boundary_conditions = boundary_condition_periodic # create homogeneous mesh -coordinates_min = -1.0 -coordinates_max = 1.0 +coordinates_min = -35.0 +coordinates_max = 35.0 N = 512 mesh = Mesh1D(coordinates_min, coordinates_max, 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 11fe5bf2..a55e72ee 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 @@ -6,12 +6,12 @@ using SparseArrays: sparse ############################################################################### # Semidiscretization of the BBM-BBM equations -equations = BBMBBMVariableEquations1D(gravity_constant = 9.81, eta0 = 2.0) +equations = BBMBBMVariableEquations1D(gravity_constant = 1.0, eta0 = 2.0) # Setup a truly discontinuous bottom topography function for this academic # testcase of well-balancedness. The errors from the analysis callback are # not important but the error for this lake-at-rest test case -# `∑|eta0-eta|` should be around machine roundoff. +# `∫|η-η₀|` should be around machine roundoff. function initial_condition_discontinuous_well_balancedness(x, t, equations::BBMBBMVariableEquations1D, mesh) @@ -38,7 +38,8 @@ N = 512 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) +accuracy_order = 4 +D1 = periodic_derivative_operator(1, accuracy_order, mesh.xmin, mesh.xmax, mesh.N) D2 = sparse(D1)^2 solver = Solver(D1, D2) diff --git a/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans.jl b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans.jl index c09c5478..b0590626 100644 --- a/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans.jl +++ b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans.jl @@ -17,7 +17,8 @@ N = 512 mesh = Mesh1D(coordinates_min, coordinates_max, N) # create solver with periodic SBP operators of accuracy order 4 -solver = Solver(mesh, 4) +accuracy_order = 4 +solver = Solver(mesh, accuracy_order) # semidiscretization holds all the necessary data structures for the spatial discretization semi = Semidiscretization(mesh, equations, initial_condition, solver, @@ -30,7 +31,7 @@ ode = semidiscretize(semi, tspan) analysis_callback = AnalysisCallback(semi; interval = 10, extra_analysis_errors = (:conservation_error,), extra_analysis_integrals = (waterheight_total, - momentum, entropy, + entropy, entropy_modified)) callbacks = CallbackSet(analysis_callback) 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 52d7e038..d84e4a45 100644 --- a/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans_relaxation.jl +++ b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans_relaxation.jl @@ -17,7 +17,8 @@ N = 512 mesh = Mesh1D(coordinates_min, coordinates_max, N) # create solver with periodic SBP operators of accuracy order 4 -solver = Solver(mesh, 4) +accuracy_order = 4 +solver = Solver(mesh, accuracy_order) # semidiscretization holds all the necessary data structures for the spatial discretization semi = Semidiscretization(mesh, equations, initial_condition, solver, 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 3f3366a1..64e0f4ad 100644 --- a/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans_upwind.jl +++ b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans_upwind.jl @@ -36,7 +36,7 @@ ode = semidiscretize(semi, tspan) analysis_callback = AnalysisCallback(semi; interval = 10, extra_analysis_errors = (:conservation_error,), extra_analysis_integrals = (waterheight_total, - momentum, entropy, + entropy, entropy_modified)) callbacks = CallbackSet(analysis_callback) 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 2844c82a..ce91e72d 100644 --- a/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_well_balanced.jl +++ b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_well_balanced.jl @@ -12,7 +12,7 @@ equations = SvaerdKalischEquations1D(gravity_constant = 1.0, eta0 = 2.0, # Setup a truly discontinuous bottom topography function for this academic # testcase of well-balancedness. The errors from the analysis callback are # not important but the error for this lake-at-rest test case -# `∑|eta0-eta|` should be around machine roundoff. +# `∫|η-η₀|` should be around machine roundoff. function initial_condition_discontinuous_well_balancedness(x, t, equations::SvaerdKalischEquations1D, mesh) @@ -39,7 +39,8 @@ N = 200 mesh = Mesh1D(coordinates_min, coordinates_max, N) # create solver with periodic SBP operators of accuracy order 4 -solver = Solver(mesh, 4) +accuracy_order = 4 +solver = Solver(mesh, accuracy_order) # semidiscretization holds all the necessary data structures for the spatial discretization semi = Semidiscretization(mesh, equations, initial_condition, solver, diff --git a/src/DispersiveShallowWater.jl b/src/DispersiveShallowWater.jl index 5757bb2b..d378b2b0 100644 --- a/src/DispersiveShallowWater.jl +++ b/src/DispersiveShallowWater.jl @@ -16,6 +16,7 @@ using SimpleUnPack: @unpack using SparseArrays: sparse using SummationByPartsOperators: AbstractDerivativeOperator, PeriodicDerivativeOperator, PeriodicUpwindOperators, + UniformPeriodicCoupledOperator, periodic_derivative_operator, derivative_order, integrate import SummationByPartsOperators: grid, xmin, xmax @@ -46,8 +47,7 @@ export Semidiscretization, semidiscretize, grid export boundary_condition_periodic -export initial_condition_convergence_test, initial_condition_sin_bathymetry, - initial_condition_dingemans +export initial_condition_convergence_test, initial_condition_dingemans export AnalysisCallback, RelaxationCallback export tstops, errors, integrals diff --git a/src/equations/bbm_bbm_1d.jl b/src/equations/bbm_bbm_1d.jl index 30446360..8e8046b5 100644 --- a/src/equations/bbm_bbm_1d.jl +++ b/src/equations/bbm_bbm_1d.jl @@ -30,7 +30,6 @@ end varnames(::typeof(prim2prim), ::BBMBBMEquations1D) = ("η", "v") varnames(::typeof(prim2cons), ::BBMBBMEquations1D) = ("h", "hv") -# TODO: Initial condition should not get a `mesh` """ initial_condition_convergence_test(x, t, equations::BBMBBMEquations1D, mesh) @@ -60,7 +59,15 @@ function create_cache(mesh, initial_condition, RealT, uEltype) - invImD2_D = (I - 1 / 6 * equations.D^2 * sparse(solver.D2)) \ Matrix(solver.D1) + if solver.D1 isa PeriodicDerivativeOperator || + solver.D1 isa UniformPeriodicCoupledOperator + invImD2_D = (I - 1 / 6 * equations.D^2 * sparse(solver.D2)) \ Matrix(solver.D1) + elseif solver.D1 isa PeriodicUpwindOperators + invImD2_D = (I - 1 / 6 * equations.D^2 * sparse(solver.D2)) \ + Matrix(solver.D1.central) + else + @error "unknown type of first-derivative operator" + end tmp1 = Array{RealT}(undef, nnodes(mesh)) return (invImD2_D = invImD2_D, tmp1 = tmp1) end diff --git a/src/equations/bbm_bbm_variable_bathymetry_1d.jl b/src/equations/bbm_bbm_variable_bathymetry_1d.jl index 784b9c9d..d69c10bf 100644 --- a/src/equations/bbm_bbm_variable_bathymetry_1d.jl +++ b/src/equations/bbm_bbm_variable_bathymetry_1d.jl @@ -30,7 +30,6 @@ end varnames(::typeof(prim2prim), ::BBMBBMVariableEquations1D) = ("η", "v", "D") varnames(::typeof(prim2cons), ::BBMBBMVariableEquations1D) = ("h", "hv", "b") -# TODO: Initial condition should not get a `mesh` """ initial_condition_convergence_test(x, t, equations::BBMBBMVariableEquations1D, mesh) @@ -59,21 +58,6 @@ function initial_condition_convergence_test(x, return SVector(eta, v, D) end -# TODO: Initial condition should not get a `mesh` -""" - initial_condition_sin_bathymetry(x, t, equations::BBMBBMVariableEquations1D, mesh) - -An initial condition with a Gaussian bump as initial water height with still water and -a sine-shaped bathymetry. -""" -function initial_condition_sin_bathymetry(x, t, equations::BBMBBMVariableEquations1D, mesh) - eta = 2.0 + 2.0 * exp(-12.0 * x^2) - v = 0.0 - D = -1.0 + 0.1 * sinpi(2.0 * x) - return SVector(eta, v, D) -end - -# TODO: Initial condition should not get a `mesh` """ initial_condition_dingemans(x, t, equations::BBMBBMVariableEquations1D, mesh) @@ -90,16 +74,16 @@ References: [link](https://repository.tudelft.nl/islandora/object/uuid:c2091d53-f455-48af-a84b-ac86680455e9/datastream/OBJ/download) """ function initial_condition_dingemans(x, t, equations::BBMBBMVariableEquations1D, mesh) - h0 = 0.8 + eta0 = 0.8 A = 0.02 - # omega = 2*pi/(2.02*sqrt(2)) - K = 0.8406220896381442 # precomputed result of find_zero(K -> omega^2 - equations.gravity * K * tanh(K * h0), 1.0) using Roots.jl - if x < -30.5 * pi / K || x > -8.5 * pi / K + # omega = 2*pi/(2.02*sqrt(2)) + k = 0.8406220896381442 # precomputed result of find_zero(k -> omega^2 - equations.gravity * k * tanh(k * eta0), 1.0) using Roots.jl + if x < -30.5 * pi / k || x > -8.5 * pi / k h = 0.0 else - h = A * cos(K * x) + h = A * cos(k * x) end - v = sqrt(equations.gravity / K * tanh(K) * h0) * h / h0 + v = sqrt(equations.gravity / k * tanh(k * eta0)) * h / eta0 if x < 11.01 || x >= 33.07 b = 0.0 elseif 11.01 <= x && x < 23.04 @@ -111,7 +95,7 @@ function initial_condition_dingemans(x, t, equations::BBMBBMVariableEquations1D, else error("should not happen") end - eta = h + h0 + eta = h + eta0 D = -b return SVector(eta, v, D) end @@ -129,7 +113,8 @@ function create_cache(mesh, D[i] = initial_condition(x[i], 0.0, equations, mesh)[3] end K = Diagonal(D .^ 2) - if solver.D1 isa PeriodicDerivativeOperator + if solver.D1 isa PeriodicDerivativeOperator || + solver.D1 isa UniformPeriodicCoupledOperator invImDKD_D = (I - 1 / 6 * sparse(solver.D1) * K * sparse(solver.D1)) \ Matrix(solver.D1) invImD2K_D = (I - 1 / 6 * sparse(solver.D2) * K) \ Matrix(solver.D1) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 57f7927d..ba3d02b3 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -81,7 +81,7 @@ of the correct length `nvariables(equations)`. """ function waterheight_total end -varnames(::typeof(waterheight_total), equations) = ("eta",) +varnames(::typeof(waterheight_total), equations) = ("η",) """ waterheight(q, equations) diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl index 68e6f945..3a53f376 100644 --- a/src/equations/svaerd_kalisch_1d.jl +++ b/src/equations/svaerd_kalisch_1d.jl @@ -44,7 +44,6 @@ end varnames(::typeof(prim2prim), ::SvaerdKalischEquations1D) = ("η", "v", "D") varnames(::typeof(prim2cons), ::SvaerdKalischEquations1D) = ("h", "hv", "b") -# TODO: Initial condition should not get a `mesh` """ initial_condition_dingemans(x, t, equations::SvaerdKalischEquations1D, mesh) @@ -64,13 +63,13 @@ function initial_condition_dingemans(x, t, equations::SvaerdKalischEquations1D, eta0 = 0.8 A = 0.02 # omega = 2*pi/(2.02*sqrt(2)) - K = 0.8406220896381442 # precomputed result of find_zero(K -> omega^2 - equations.gravity * K * tanh(K * eta0), 1.0) using Roots.jl - if x < -30.5 * pi / K || x > -8.5 * pi / K + k = 0.8406220896381442 # precomputed result of find_zero(k -> omega^2 - equations.gravity * k * tanh(k * eta0), 1.0) using Roots.jl + if x < -30.5 * pi / k || x > -8.5 * pi / k h = 0.0 else - h = A * cos(K * x) + h = A * cos(k * x) end - v = sqrt(equations.gravity / K * tanh(K) * eta0) * h / eta0 + v = sqrt(equations.gravity / k * tanh(k * eta0)) * h / eta0 if x < 11.01 || x >= 33.07 b = 0.0 elseif 11.01 <= x && x < 23.04 @@ -108,7 +107,8 @@ function create_cache(mesh, tmp1 = similar(h) tmp2 = similar(h) hmD1betaD1 = Array{RealT}(undef, nnodes(mesh), nnodes(mesh)) - if solver.D1 isa PeriodicDerivativeOperator + if solver.D1 isa PeriodicDerivativeOperator || + solver.D1 isa UniformPeriodicCoupledOperator D1_central = solver.D1 sparse_D1 = sparse(D1_central) D1betaD1 = sparse_D1 * Diagonal(beta_hat) * sparse_D1 @@ -142,7 +142,8 @@ function rhs!(du_ode, u_ode, t, mesh, equations::SvaerdKalischEquations1D, @. h = eta + D hv = h .* v - if solver.D1 isa PeriodicDerivativeOperator + if solver.D1 isa PeriodicDerivativeOperator || + solver.D1 isa UniformPeriodicCoupledOperator D1eta = D1_central * eta D1v = D1_central * v tmp1 = alpha_hat .* (D1_central * (alpha_hat .* D1eta)) @@ -170,6 +171,12 @@ function rhs!(du_ode, u_ode, t, mesh, equations::SvaerdKalischEquations1D, 0.5 * (vD1y - D1vy - yD1v) - 0.5 * D1_central * (gamma_hat .* (solver.D2 * v)) - 0.5 * solver.D2 * (gamma_hat .* D1v)) + # not split form + # tmp2 = -(D1_central * (hv .* v) - v .* (D1_central * hv)+ + # equations.gravity * h .* D1eta + + # vD1y - D1vy - + # 0.5 * D1_central * (gamma_hat .* (solver.D2 * v)) - + # 0.5 * solver.D2 * (gamma_hat .* D1v)) dv[:] = hmD1betaD1 \ tmp2 return nothing @@ -238,7 +245,8 @@ number of nodes as length of the second dimension. v = view(q, 2, :) D = view(q, 3, :) beta_hat = equations.beta * (eta .+ D) .^ 3 - if cache.D1 isa PeriodicDerivativeOperator + if cache.D1 isa PeriodicDerivativeOperator || + cache.D1 isa UniformPeriodicCoupledOperator tmp = 0.5 * beta_hat .* ((cache.D1 * v) .^ 2) elseif cache.D1 isa PeriodicUpwindOperators tmp = 0.5 * beta_hat .* ((cache.D1.minus * v) .^ 2) diff --git a/src/util.jl b/src/util.jl index c4267ac5..df7d7378 100644 --- a/src/util.jl +++ b/src/util.jl @@ -61,7 +61,7 @@ include(DispersiveShallowWater.path_create_figures()) ``` """ function path_create_figures() - pkgdir(DispersiveShallowWater, "create_figures.jl") + pkgdir(DispersiveShallowWater, "visualization", "create_figures.jl") end # Note: We can't call the method below `DispersiveShallowWater.include` since that is created automatically diff --git a/src/visualization.jl b/src/visualization.jl index 62567773..f80ade23 100644 --- a/src/visualization.jl +++ b/src/visualization.jl @@ -52,7 +52,7 @@ end end plot_title --> "$(get_name(semi.equations)) at t = $(round(t, digits=5))" - layout := nsubplots + layout --> nsubplots for i in 1:nsubplots # Don't plot bathymetry in separate subplot @@ -60,7 +60,7 @@ end if plot_initial == true @series begin - subplot := i + subplot --> i linestyle := :solid label := "initial $(names[i])" grid(semi), q_exact[i, :] @@ -68,7 +68,7 @@ end end @series begin - subplot := i + subplot --> i label --> names[i] xguide --> "x" yguide --> names[i] @@ -80,12 +80,12 @@ end # Plot the bathymetry if plot_bathymetry == true @series begin - subplot := 1 + subplot --> 1 linestyle := :solid label := "bathymetry" - xguide := "x" - yguide := names[1] - title := names[1] + xguide --> "x" + yguide --> names[1] + title --> names[1] color := :black grid(semi), bathy end @@ -121,18 +121,17 @@ end names = varnames(conversion, equations) plot_title --> "$(get_name(semi.equations)) at x = $(round(grid(semi)[index], digits=5))" - size --> (1200, 800) - layout := nsubplots + layout --> nsubplots for i in 1:nsubplots # Don't plot bathymetry in separate subplot names[i] in ["D", "b"] && continue @series begin - subplot := i - label := names[i] - xguide := "t" - yguide := names[i] - title := names[i] + subplot --> i + label --> names[i] + xguide --> "t" + yguide --> names[i] + title --> names[i] sol.t, data[i, :] end end @@ -167,12 +166,16 @@ function pretty(name) return "∫P" elseif name == :entropy return "∫U" + elseif name == :entropy_modified + return "∫U_mod" elseif name == :l2_error return "L² error" elseif name == :linf_error return "L∞ error" elseif name == :conservation_error return "∫|q_q₀|" + elseif name == :lake_at_rest_error + return "∫|η-η₀|" else return string(name) end @@ -206,7 +209,7 @@ end name in exclude && continue @series begin subplot --> subplot - label := pretty(name) * " " * label_extension + label --> pretty(name) * " " * label_extension title --> "errors" xguide --> "t" yguide --> "sum of errors" diff --git a/test/test_bbm_bbm_variable_bathymetry_1d.jl b/test/test_bbm_bbm_variable_bathymetry_1d.jl index 8aca8dfe..8bcea168 100644 --- a/test/test_bbm_bbm_variable_bathymetry_1d.jl +++ b/test/test_bbm_bbm_variable_bathymetry_1d.jl @@ -25,8 +25,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_variable_bathymetry_1d") @test_trixi_include(joinpath(EXAMPLES_DIR, "bbm_bbm_variable_bathymetry_1d_relaxation.jl"), tspan=(0.0, 1.0), - l2=[1.234497154595864 2.0325352963894225 0.0], - linf=[1.8566998050289607 2.537410906888876 0.0], + l2=[0.001695855504730616 0.0034287779832081213 0.0], + linf=[0.0014540429957374812 0.001806888012062302 0.0], cons_error=[2.0694557179012918e-13 2.789919043027633e-13 0.0], change_waterheight=-2.1271873151818e-13, change_velocity=2.7874259106214573e-13, @@ -37,8 +37,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_variable_bathymetry_1d") @test_trixi_include(joinpath(EXAMPLES_DIR, "bbm_bbm_variable_bathymetry_1d_dg_upwind_relaxation.jl"), tspan=(0.0, 1.0), - l2=[1.2344971565355631 2.0325352892395023 0.0], - linf=[1.8563719900575752 2.537410881290326 0.0], + l2=[0.004960390675877118 0.00830158298558452 0.0], + linf=[0.011473564092181476 0.016650286103717254 0.0], cons_error=[3.1796787425264483e-13 1.080125572316959e-14 0.0], change_waterheight=3.1796787425264483e-13, change_velocity=-1.0584519372081047e-12, @@ -49,8 +49,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_variable_bathymetry_1d") @test_trixi_include(joinpath(EXAMPLES_DIR, "bbm_bbm_variable_bathymetry_1d_upwind_relaxation.jl"), tspan=(0.0, 1.0), - l2=[1.2344971597115844 2.032535283544955 0.0], - linf=[1.8566997938473593 2.5374108392574604 0.0], + l2=[0.002644536848587577 0.005179311522968244 0.0], + linf=[0.0022397723587049834 0.0031376618593412786 0.0], cons_error=[7.194245199571014e-14 1.2388735870505485e-12 0.0], change_waterheight=-1.5765166949677223e-13, change_velocity=1.8791999206735355e-13, @@ -75,12 +75,12 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_variable_bathymetry_1d") "bbm_bbm_variable_bathymetry_1d_dingemans.jl"), tspan=(0.0, 1.0), N=512, - l2=[0.23119432844001495 0.7543986688031077 0.0], - linf=[0.037067680051741214 0.12128735996476576 0.0], + l2=[0.2322073114427607 0.7753458687737584 0.0], + linf=[0.037222719015511885 0.124336213226626 0.0], cons_error=[1.4210854715202004e-13 3.1478183774857893e-15 0.0], change_waterheight=-1.4210854715202004e-13, change_velocity=-3.1478183774857893e-15, - change_entropy=-1.3967564882477745e-9) + change_entropy=-1.442231223336421e-9) end end diff --git a/test/test_svaerd_kalisch_1d.jl b/test/test_svaerd_kalisch_1d.jl index 793829f3..aec84ee6 100644 --- a/test/test_svaerd_kalisch_1d.jl +++ b/test/test_svaerd_kalisch_1d.jl @@ -13,26 +13,24 @@ EXAMPLES_DIR = joinpath(examples_dir(), "svaerd_kalisch_1d") "svaerd_kalisch_1d_dingemans.jl"), tspan=(0.0, 1.0), N=512, - l2=[0.22648616606810756 0.7325813437573053 0.0], - linf=[0.03645879748792702 0.11849678361651249 0.0], - cons_error=[3.979039320256561e-13 7.670877062280329e-5 0.0], + l2=[0.22796106962338855 0.7519327063662515 0.0], + linf=[0.036708347831218346 0.12141172207472928 0.0], + cons_error=[3.979039320256561e-13 4.937137540373564e-5 0.0], change_waterheight=-3.979039320256561e-13, - change_momentum=-2.109259079896564e-9, - change_entropy=-0.0005910026394531087, - change_entropy_modified=3.20516357987799e-6) + change_entropy=-0.00024362648639453255, + change_entropy_modified=3.240868750253867e-6) end @trixi_testset "svaerd_kalisch_1d_dingemans_upwind" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "svaerd_kalisch_1d_dingemans_upwind.jl"), tspan=(0.0, 1.0), - l2=[0.22656501380746077 0.7327645752975466 0.0], - linf=[0.03648151661791921 0.1182504804050214 0.0], - cons_error=[1.1368683772161603e-13 7.59953335162671e-5 0.0], + l2=[0.2280370308166863 0.7521344942401095 0.0], + linf=[0.03673101553812019 0.12116306036094074 0.0], + cons_error=[1.1368683772161603e-13 4.871598417571836e-5 0.0], change_waterheight=-1.1368683772161603e-13, - change_momentum=-2.28489619585881e-9, - change_entropy=-0.0005840425510541536, - change_entropy_modified=2.108449052684591e-6) + change_entropy=-0.00023645232727176335, + change_entropy_modified=2.0682244894487667e-6) end @trixi_testset "svaerd_kalisch_1d_dingemans_relaxation" begin @@ -40,11 +38,11 @@ EXAMPLES_DIR = joinpath(examples_dir(), "svaerd_kalisch_1d") "svaerd_kalisch_1d_dingemans_relaxation.jl"), tspan=(0.0, 1.0), N=512, - l2=[0.22648560566941997 0.7325790781930053 0.0], - linf=[0.03645867682442583 0.1184963632507331 0.0], - cons_error=[3.979039320256561e-13 7.670526860626775e-5 0.0], + l2=[0.2279601359992612 0.7519295302772953 0.0], + linf=[0.03670816991272552 0.12141112705223758 0.0], + cons_error=[3.979039320256561e-13 4.937017087502672e-5 0.0], change_waterheight=-3.979039320256561e-13, - change_entropy=-0.0005940058100577517, + change_entropy=-0.0002466715038735856, change_entropy_modified=0.0) end diff --git a/test/test_unit.jl b/test/test_unit.jl index 3b8062d9..80b0acaf 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -33,11 +33,11 @@ using SparseArrays: sparse, SparseMatrixCSC @test real(solver) == Float64 @test_nowarn show(stdout, solver) - Dop = legendre_derivative_operator(-1.0, 1.0, p + 1) - sbp_mesh = UniformPeriodicMesh1D(-1.0, 1.0, 512 ÷ (p + 1)) - central = couple_discontinuously(Dop, sbp_mesh) - minus = couple_discontinuously(Dop, sbp_mesh, Val(:minus)) - plus = couple_discontinuously(Dop, sbp_mesh, Val(:plus)) + D_legendre = legendre_derivative_operator(-1.0, 1.0, p + 1) + uniform_mesh = UniformPeriodicMesh1D(-1.0, 1.0, 512 ÷ (p + 1)) + central = couple_discontinuously(D_legendre, uniform_mesh) + minus = couple_discontinuously(D_legendre, uniform_mesh, Val(:minus)) + plus = couple_discontinuously(D_legendre, uniform_mesh, Val(:plus)) D2 = sparse(plus) * sparse(minus) solver = @test_nowarn Solver(central, D2) @test solver.D1 isa UniformPeriodicCoupledOperator diff --git a/visualization/bbm_bbm_1d_widestencil.jl b/visualization/bbm_bbm_1d_widestencil.jl new file mode 100644 index 00000000..8dd057b0 --- /dev/null +++ b/visualization/bbm_bbm_1d_widestencil.jl @@ -0,0 +1,43 @@ +using OrdinaryDiffEq +using DispersiveShallowWater +using SummationByPartsOperators: periodic_derivative_operator +using SparseArrays: sparse + +############################################################################### +# Semidiscretization of the BBM-BBM equations + +equations = BBMBBMEquations1D(gravity_constant = 9.81, D = 2.0) + +# initial_condition_convergence_test needs periodic boundary conditions +initial_condition = initial_condition_convergence_test +boundary_conditions = boundary_condition_periodic + +# create homogeneous mesh +coordinates_min = -35.0 +coordinates_max = 35.0 +N = 512 +mesh = Mesh1D(coordinates_min, coordinates_max, N) + +# create solver with periodic SBP operators of accuracy order 4 +accuracy_order = 4 +D1 = periodic_derivative_operator(1, accuracy_order, mesh.xmin, mesh.xmax, mesh.N) +D2 = sparse(D1)^2 +solver = Solver(D1, D2) + +# semidiscretization holds all the necessary data structures for the spatial discretization +semi = Semidiscretization(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# Create `ODEProblem` and run the simulation +tspan = (0.0, 30.0) +ode = semidiscretize(semi, tspan) +analysis_callback = AnalysisCallback(semi; interval = 10, + extra_analysis_errors = (:conservation_error,), + extra_analysis_integrals = (waterheight_total, + velocity, entropy)) +callbacks = CallbackSet(analysis_callback) + +saveat = range(tspan..., length = 100) +sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7, + save_everystep = false, callback = callbacks, saveat = saveat) diff --git a/visualization/bbm_bbm_variable_bathymetry_1d_upwind.jl b/visualization/bbm_bbm_variable_bathymetry_1d_upwind.jl new file mode 100644 index 00000000..3ab09684 --- /dev/null +++ b/visualization/bbm_bbm_variable_bathymetry_1d_upwind.jl @@ -0,0 +1,45 @@ +using OrdinaryDiffEq +using DispersiveShallowWater +using SummationByPartsOperators: upwind_operators, periodic_derivative_operator +using SparseArrays: sparse + +############################################################################### +# Semidiscretization of the BBM-BBM equations + +equations = BBMBBMVariableEquations1D(gravity_constant = 9.81) + +# initial_condition_convergence_test needs periodic boundary conditions +initial_condition = initial_condition_convergence_test +boundary_conditions = boundary_condition_periodic + +# create homogeneous mesh +coordinates_min = -35.0 +coordinates_max = 35.0 +N = 512 +mesh = Mesh1D(coordinates_min, coordinates_max, N) + +# create solver with periodic SBP operators of accuracy order 4 +accuracy_order = 4 +D1 = upwind_operators(periodic_derivative_operator; derivative_order = 1, + accuracy_order = accuracy_order, xmin = mesh.xmin, xmax = mesh.xmax, + N = mesh.N) +D2 = sparse(D1.plus) * sparse(D1.minus) +solver = Solver(D1, D2) + +# semidiscretization holds all the necessary data structures for the spatial discretization +semi = Semidiscretization(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# Create `ODEProblem` and run the simulation +tspan = (0.0, 30.0) +ode = semidiscretize(semi, tspan) +analysis_callback = AnalysisCallback(semi; interval = 10, + extra_analysis_errors = (:conservation_error,), + extra_analysis_integrals = (waterheight_total, + velocity, entropy)) +callbacks = CallbackSet(analysis_callback) + +saveat = range(tspan..., length = 100) +sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7, + save_everystep = false, callback = callbacks, saveat = saveat) diff --git a/visualization/bbm_bbm_variable_bathymetry_1d_widestencil.jl b/visualization/bbm_bbm_variable_bathymetry_1d_widestencil.jl new file mode 100644 index 00000000..8f62c419 --- /dev/null +++ b/visualization/bbm_bbm_variable_bathymetry_1d_widestencil.jl @@ -0,0 +1,43 @@ +using OrdinaryDiffEq +using DispersiveShallowWater +using SummationByPartsOperators: periodic_derivative_operator +using SparseArrays: sparse + +############################################################################### +# Semidiscretization of the BBM-BBM equations + +equations = BBMBBMVariableEquations1D(gravity_constant = 9.81) + +# initial_condition_convergence_test needs periodic boundary conditions +initial_condition = initial_condition_convergence_test +boundary_conditions = boundary_condition_periodic + +# create homogeneous mesh +coordinates_min = -35.0 +coordinates_max = 35.0 +N = 512 +mesh = Mesh1D(coordinates_min, coordinates_max, N) + +# create solver with periodic SBP operators of accuracy order 4 +accuracy_order = 4 +D1 = periodic_derivative_operator(1, accuracy_order, mesh.xmin, mesh.xmax, mesh.N) +D2 = sparse(D1)^2 +solver = Solver(D1, D2) + +# semidiscretization holds all the necessary data structures for the spatial discretization +semi = Semidiscretization(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# Create `ODEProblem` and run the simulation +tspan = (0.0, 30.0) +ode = semidiscretize(semi, tspan) +analysis_callback = AnalysisCallback(semi; interval = 10, + extra_analysis_errors = (:conservation_error,), + extra_analysis_integrals = (waterheight_total, + velocity, entropy)) +callbacks = CallbackSet(analysis_callback) + +saveat = range(tspan..., length = 100) +sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7, + save_everystep = false, callback = callbacks, saveat = saveat) diff --git a/visualization/create_figures.jl b/visualization/create_figures.jl new file mode 100644 index 00000000..b797d337 --- /dev/null +++ b/visualization/create_figures.jl @@ -0,0 +1,733 @@ +using DispersiveShallowWater +using SummationByPartsOperators: SummationByPartsOperators, + periodic_derivative_operator, + upwind_operators, + legendre_derivative_operator, + legendre_second_derivative_operator, + UniformPeriodicMesh1D, + couple_discontinuously, + couple_continuously +using SparseArrays: sparse +using Plots +using LaTeXStrings + +const OUT = "out/" +ispath(OUT) || mkpath(OUT) +const VISUALIZATION_DIR = pkgdir(DispersiveShallowWater, "visualization") +const EXAMPLES_DIR_BBMBBM = joinpath(examples_dir(), "bbm_bbm_1d") +const EXAMPLES_DIR_BBMBBM_VARIABLE = joinpath(examples_dir(), + "bbm_bbm_variable_bathymetry_1d") +const EXAMPLES_DIR_SVAERD_KALISCH = joinpath(examples_dir(), "svaerd_kalisch_1d") + +# Chapter 2 +# Plot of bathymetry and waterheight +function fig_1() + L = 1.0 + n = 100 + x = LinRange(0.0, L, n) + fontsize = 20 + + # just pick some function for b and eta that look nice + H = 1.012 + + b(x) = x * cos.(3 * pi * x) + H + plot(x, b, color = :gray, fill = (0, 0.8, :gray), fillstyle = :/, linewidth = 3, + legend = nothing, ticks = nothing, border = :none) + + eta(x) = x / (x^2 + 1) * sin(2 * pi * x) + H + 1.5 + plot!(x, eta, color = :blue, fill = (b.(x), 0.4, :blue), linewidth = 3) + + x1 = 0.2 + plot!([x1, x1], [b(x1), eta(x1)], line = (Plots.Arrow(:open, :both, 2.5, 2.0), :black), + annotation = (x1 - 0.08, (eta(x1) + b(x1)) / 2, text(L"h(t, x)", fontsize)), + linewidth = 2) + x2 = 0.4 + plot!([x2, x2], [0.0, b(x2)], line = (Plots.Arrow(:open, :both, 2.5, 2.0), :black), + annotation = (x2 + 0.06, b(x2) / 2, text(L"b(x)", fontsize)), linewidth = 2) + x3 = 0.8 + plot!([x3, x3], [0.0, eta(x3)], line = (Plots.Arrow(:open, :both, 2.5, 2.0), :black), + annotation = (x3 - 0.08, eta(x3) / 2, text(L"\eta(t, x)", fontsize)), + linewidth = 2) + + savefig(joinpath(OUT, "bathymetry.pdf")) +end + +# Plot of dispersion relations +function fig_2() + linewidth = 2 + markersize = 5 + + h0 = 1.0 + g = 1.0 + c0 = sqrt(g * h0) + + k = 0.01:0.5:(8 * pi) + k_zoom = 0.01:0.3:pi + ylim = (0.0, 1.1) + + omega_euler(k) = sqrt(g * k) * sqrt(tanh(h0 * k)) + c_euler(k) = omega_euler(k) / k + 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, 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, + framestyle = :box) + + function plot_dispersion_relation(omega, label, markershape) + c(k) = omega(k) / k + plot!(k, c.(k) ./ c0, label = label, linewidth = linewidth, + markershape = markershape, markersize = markersize) + plot!(k_zoom, c.(k_zoom) ./ c0, subplot = 2, linewidth = linewidth, + markershape = markershape, markersize = markersize) + end + + omega_bbmbbm_(k, d0) = sqrt(g * h0) * k / (1 + 1 / 6 * (d0 * k)^2) + omega_bbmbbm(k) = omega_bbmbbm_(k, h0) + plot_dispersion_relation(omega_bbmbbm, "BBM-BBM", :cross) + + alpha_set1 = -1 / 3 * c0 * h0^2 + beta_set1 = 0.0 * h0^3 + gamma_set1 = 0.0 * c0 * h0^3 + + alpha_set2 = 0.0004040404040404049 * c0 * h0^2 + beta_set2 = 0.49292929292929294 * h0^3 + gamma_set2 = 0.15707070707070708 * c0 * h0^3 + + alpha_set3 = 0.0 * c0 * h0^2 + beta_set3 = 0.27946992481203003 * h0^3 + gamma_set3 = 0.0521077694235589 * c0 * h0^3 + + alpha_set4 = 0.0 * c0 * h0^2 + beta_set4 = 0.2308939393939394 * h0^3 + gamma_set4 = 0.04034343434343434 * c0 * h0^3 + + function char_equation(alpha, beta, gamma, k) + a = (1 + beta / h0 * k^2) + b = (-alpha - beta * alpha / h0 * k^2 - gamma / h0) * k^3 + c = -g * h0 * k^2 + gamma * alpha / h0 * k^6 + omega1 = (-b + sqrt(b^2 - 4 * a * c)) / (2 * a) + # omega2 = (-b - sqrt(b^2 - 4*a*c))/(2*a) + return omega1 + end + + omega_set1(k) = char_equation(alpha_set1, beta_set1, gamma_set1, k) + plot_dispersion_relation(omega_set1, "S.-K. set 1", :rtriangle) + + omega_set2(k) = char_equation(alpha_set2, beta_set2, gamma_set2, k) + plot_dispersion_relation(omega_set2, "S.-K. set 2", :star5) + + omega_set3(k) = char_equation(alpha_set3, beta_set3, gamma_set3, k) + plot_dispersion_relation(omega_set3, "S.-K. set 3", :star8) + + omega_set4(k) = char_equation(alpha_set4, beta_set4, gamma_set4, k) + plot_dispersion_relation(omega_set4, "S.-K. set 4", :diamond) + + # Plot box + plot!([0.0, pi], [0.54, 0.54], color = :black, label = :none) + plot!([0.0, pi], [1.0, 1.0], color = :black, label = :none) + plot!([0.0, 0.0], [0.54, 1.0], color = :black, label = :none) + plot!([pi, pi], [0.54, 1.0], color = :black, label = :none) + + # Plot connecting lines + plot!([pi, 6.8], [0.54, 0.629], color = :black, label = :none) + plot!([pi, 6.8], [1, 1.01], color = :black, label = :none) + + savefig(joinpath(OUT, "dispersion_relations.pdf")) +end + +# Chapter 4 +# Chapter 4.1 Soliton + +const OUT_SOLITON = joinpath(OUT, "soliton") +ispath(OUT_SOLITON) || mkpath(OUT_SOLITON) + +# Plot convergence orders for baseline and relaxation +function fig_3() + tspan = (0.0, 10.0) + accuracy_orders = [2, 4, 6, 8] + iters = [4, 4, 4, 3] + initial_Ns = [128, 128, 128, 128] + + all_Ns = minimum(initial_Ns) * 2 .^ (0:(maximum(iters) - 1)) + + 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(joinpath(EXAMPLES_DIR_BBMBBM, + "bbm_bbm_1d_basic.jl"), + iters[i]; N = initial_Ns[i], tspan = tspan, + accuracy_order = accuracy_orders[i]) + # Use sum over all L^2-errors for all variables, i.e. ||η - η_ana||_2 + ||v - v_ana||_2 + 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, label = "p = $(accuracy_orders[i]), EOC: $eoc_mean", + markershape = markershapes[i], linewidth = linewidth, markersize = markersize, + subplot = 1) + end + + # right subplot: relaxation + for i in 1:length(accuracy_orders) + Ns = initial_Ns[i] * 2 .^ (0:(iters[i] - 1)) + _, errormatrix = convergence_test(joinpath(EXAMPLES_DIR_BBMBBM, + "bbm_bbm_1d_relaxation.jl"), + iters[i]; N = initial_Ns[i], tspan = tspan, + accuracy_order = accuracy_orders[i]) + # Use sum over all L^2-errors for all variables, i.e. ||η - η_ana||_2 + ||v - v_ana||_2 + 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, label = "p = $(accuracy_orders[i]), EOC: $eoc_mean", + markershape = markershapes[i], linewidth = linewidth, markersize = markersize, + subplot = 2) + end + xticks!(all_Ns, string.(all_Ns)) + savefig(joinpath(OUT_SOLITON, "orders.pdf")) +end + +# Plot errors, change of invariants, and solution at final time for baseline and relaxation +function fig_4_5_6() + linewidth = 2 + linestyles = [:dash, :dot] + + g = 9.81 + D = 2.0 + c = 5 / 2 * sqrt(g * D) + xmin = -35.0 + xmax = 35.0 + tspan = (0.0, 50 * (xmax - xmin) / c) + N = 512 + accuracy_order = 8 + + # baseline + trixi_include(joinpath(EXAMPLES_DIR_BBMBBM, "bbm_bbm_1d_basic.jl"), + gravity_constant = g, D = D, coordinates_min = xmin, + coordinates_max = xmax, tspan = tspan, N = N, + accuracy_order = accuracy_order) + p1 = plot(analysis_callback, title = "", label_extension = "baseline", + linestyles = [:solid :dash :dot], + linewidth = linewidth, layout = 2, subplot = 1) + p2 = plot(analysis_callback, title = "", what = (:errors,), + label_extension = "baseline", linestyle = linestyles[1], + 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, plot_bathymetry = false, + linestyle = linestyles[1], linewidth = linewidth, plot_title = "", title = "", + ylims = [(-8, 3) (-1, 40)]) + x = DispersiveShallowWater.grid(semi) + q = DispersiveShallowWater.wrap_array(sol.u[end], semi) + plot!(p3, x, view(q, 1, :), inset = (1, bbox(0.11, 0.6, 0.35, 0.32)), subplot = 3, + xlim = (-20, -10), + ylim = (-0.05, 0.05), legend = nothing, linewidth = linewidth, + linestyle = linestyles[1], + color = 2, + tickfontsize = 5, yticks = [0.04, 0.0, -0.04], xticks = [-20, -15, -10], + plot_initial = true, plot_bathymetry = false, framestyle = :box) + q_exact = DispersiveShallowWater.wrap_array(DispersiveShallowWater.compute_coefficients(initial_condition, + tspan[2], + semi), + semi) + plot!(p3, x, view(q_exact, 1, :), subplot = 3, legend = nothing, linewidth = linewidth, + linestyle = :solid, color = 1) + + # relaxation + trixi_include(joinpath(EXAMPLES_DIR_BBMBBM, "bbm_bbm_1d_relaxation.jl"), + gravity_constant = g, D = D, coordinates_min = xmin, + coordinates_max = xmax, tspan = tspan, N = N, + accuracy_order = accuracy_order) + plot!(p1, analysis_callback, title = "", label_extension = "relaxation", + linestyles = [:solid :dash :dot], + linewidth = linewidth, subplot = 2) + plot!(p2, analysis_callback, title = "", what = (:errors,), + label_extension = "relaxation", linestyle = linestyles[2], 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 = linestyles[2], + linewidth = linewidth, plot_title = "", title = "", color = 3) + x = DispersiveShallowWater.grid(semi) + q = DispersiveShallowWater.wrap_array(sol.u[end], semi) + plot!(p3, x, view(q, 1, :), subplot = 3, legend = nothing, linewidth = linewidth, + linestyle = linestyles[2], color = 3) + + # Plot box + plot!(p3, [-20, -10], [-0.1, -0.1], color = :black, label = :none) + plot!(p3, [-20, -10], [0.1, 0.1], color = :black, label = :none) + plot!(p3, [-20, -20], [-0.1, 0.1], color = :black, label = :none) + plot!(p3, [-10, -10], [-0.1, 0.1], color = :black, label = :none) + + # Plot connecting lines + plot!(p3, [-20, -29], [-0.1, -3.6], color = :black, label = :none) + plot!(p3, [-10, -3.15], [-0.1, -3.6], color = :black, label = :none) + + savefig(p1, joinpath(OUT_SOLITON, "invariants.pdf")) + savefig(p2, joinpath(OUT_SOLITON, "errors.pdf")) + savefig(p3, joinpath(OUT_SOLITON, "solution.pdf")) +end + +# Plot errors for narrow-stencil, wide-stencil and upwind operators (all using relaxation) +function fig_7() + linewidth = 2 + linestyles = [:solid, :dash, :dot, :dashdot] + + g = 9.81 + D = 2.0 + c = 5 / 2 * sqrt(g * D) + xmin = -35.0 + xmax = 35.0 + tspan = (0.0, 20 * (xmax - xmin) / c) + N = 512 + accuracy_order = 8 + + plot() + + D1 = periodic_derivative_operator(1, accuracy_order, xmin, xmax, N) + D2 = sparse(D1)^2 + solver_widestencil = Solver(D1, D2) + + D1 = periodic_derivative_operator(1, accuracy_order, xmin, xmax, N) + D2 = periodic_derivative_operator(2, accuracy_order, xmin, xmax, N) + solver_narrowstencil = Solver(D1, D2) + + D1 = upwind_operators(periodic_derivative_operator; derivative_order = 1, + accuracy_order = accuracy_order, xmin = xmin, xmax = xmax, + N = N) + D2 = sparse(D1.plus) * sparse(D1.minus) + solver_upwind = Solver(D1, D2) + solvers = [ + solver_narrowstencil, + solver_narrowstencil, + solver_widestencil, + solver_upwind, + ] + labels = [ + "narrow-stencil", + "narrow-stencil in velocity equation", + "wide-stencil", + "upwind", + ] + examples = [joinpath(EXAMPLES_DIR_BBMBBM, "bbm_bbm_1d_relaxation.jl"), + joinpath(EXAMPLES_DIR_BBMBBM_VARIABLE, + "bbm_bbm_variable_bathymetry_1d_relaxation.jl"), + joinpath(EXAMPLES_DIR_BBMBBM_VARIABLE, + "bbm_bbm_variable_bathymetry_1d_relaxation.jl"), + joinpath(EXAMPLES_DIR_BBMBBM_VARIABLE, + "bbm_bbm_variable_bathymetry_1d_relaxation.jl")] + + for (i, solver) in enumerate(solvers) + trixi_include(examples[i], + gravity_constant = g, D = D, coordinates_min = xmin, + coordinates_max = xmax, tspan = tspan, N = N, + accuracy_order = accuracy_order, solver = solver) + plot!(analysis_callback, title = "", what = (:errors,), + label_extension = labels[i], linestyle = linestyles[i], + linewidth = linewidth, + ylabel = L"\Vert\eta - \eta_{ana}\Vert_2 + \Vert v - v_{ana}\Vert_2", + exclude = [:conservation_error, :linf_error]) + end + + savefig(joinpath(OUT_SOLITON, "errors_stencils.pdf")) +end + +# Compare the orders of narrow-stencil, wide-stencil and upwind SBP operators +# Not used in the thesis, but nonetheless interesting +function fig_orders_different_stencils() + tspan = (0.0, 10.0) + xmin = -35.0 + xmax = 35.0 + accuracy_orders = [2, 4, 6, 8] + iters = [4, 4, 4, 3] + initial_Ns = [128, 128, 128, 128] + + all_Ns = minimum(initial_Ns) * 2 .^ (0:(maximum(iters) - 1)) + + linewidth = 2 + markersize = 5 + markershapes = [:circle, :star5, :star8, :rtriangle] + plot(label = :none, xscale = :log2, yscale = :log10, xlabel = "N", ylim = (1e-4, 1e3), + ylabel = L"\Vert\eta - \eta_{ana}\Vert_2 + \Vert v - v_{ana}\Vert_2", + legend = :topright, layout = (1, 3)) + + # put examples in separate files since the different solvers cannot be set with the convergence_test + # because they depend on N + examples = [ + joinpath(EXAMPLES_DIR_BBMBBM_VARIABLE, "bbm_bbm_variable_bathymetry_1d_basic.jl"), + joinpath(VISUALIZATION_DIR, "bbm_bbm_variable_bathymetry_1d_widestencil.jl"), + joinpath(VISUALIZATION_DIR, "bbm_bbm_variable_bathymetry_1d_upwind.jl")] + for (j, example) in enumerate(examples) + for i in 1:length(accuracy_orders) + Ns = initial_Ns[i] * 2 .^ (0:(iters[i] - 1)) + _, errormatrix = convergence_test(example, + iters[i]; N = initial_Ns[i], tspan = tspan, + accuracy_order = accuracy_orders[i], + coordinates_min = xmin, + coordinates_max = xmax) + # Use sum over all L^2-errors for all variables, i.e. ||η - η_ana||_2 + ||v - v_ana||_2 + 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, label = "p = $(accuracy_orders[i]), EOC: $eoc_mean", + markershape = markershapes[i], linewidth = linewidth, + markersize = markersize, + subplot = j) + end + end + + xticks!(all_Ns, string.(all_Ns)) + plot!(top_margin = 3 * Plots.mm, subplot = 1) + savefig(joinpath(OUT_SOLITON, "orders_stencils.pdf")) +end + +# Chapter 4.2 Lake-at-rest +const OUT_LAKEATREST = joinpath(OUT, "lake_at_rest") +ispath(OUT_LAKEATREST) || mkpath(OUT_LAKEATREST) + +# Lake-at-rest error for long-time simulation with discontinuous bottom +function fig_8() + linewidth = 2 + N = 200 + accuracy_order = 4 + xmin = -1.0 + xmax = 1.0 + tspan = (0.0, 100.0) + D1 = upwind_operators(periodic_derivative_operator; derivative_order = 1, + accuracy_order = accuracy_order, xmin = xmin, xmax = xmax, + N = N) + D2 = sparse(D1.plus) * sparse(D1.minus) + solver = Solver(D1, D2) + trixi_include(joinpath(EXAMPLES_DIR_BBMBBM_VARIABLE, + "bbm_bbm_variable_bathymetry_1d_well_balanced.jl"); + N = N, tspan = tspan, solver = solver, dt = 0.5) + plot(analysis_callback, exclude = [:waterheight_total, :velocity, :entropy], + label_extension = "BBM-BBM", plot_title = "", title = "", + ylabel = "lake-at-rest error", linewidth = linewidth) + + trixi_include(joinpath(EXAMPLES_DIR_SVAERD_KALISCH, + "svaerd_kalisch_1d_well_balanced.jl"); + N = N, tspan = tspan, solver = solver, dt = 0.003) + plot!(analysis_callback, exclude = [:waterheight_total, :momentum, :entropy], + label_extension = "Svärd-Kalisch", plot_title = "", title = "", + ylabel = "lake-at-rest error", linewidth = linewidth) + savefig(joinpath(OUT_LAKEATREST, "lake_at_rest_error_discontinuous.pdf")) +end + +# Plot of condition number of matrix that needs to be inverted for the Svärd-Kalisch equations for different order of accuracy +# Not used in the thesis, but nonetheless interesting +function fig_condition_number() + xmin = -1.0 + xmax = 1.0 + accuracy_orders = [2, 4, 6, 8] + eta0 = 2.0 + beta = 0.49292929292929294 + Ns = 10:10:500 + conds = Array{Float64}(undef, length(Ns)) + plot(xlabel = "N", ylabel = L"cond_2") + for accuracy_order in accuracy_orders + for (i, N) in enumerate(Ns) + D1 = upwind_operators(periodic_derivative_operator; derivative_order = 1, + accuracy_order = accuracy_order, xmin = xmin, xmax = xmax, + N = N) + eta = fill(eta0, N) + D = fill(-1.0, N) + for (i, x) in enumerate(SummationByPartsOperators.grid(D1)) + if x >= 0.5 && x <= 0.75 + D[i] = -1.5 - 0.5 * sinpi(2.0 * x) + end + end + d = eta0 .+ D + beta_hat = beta * d .^ 3 + + D1betaD1 = sparse(D1.plus) * Diagonal(beta_hat) * sparse(D1.minus) + hmD1betaD1 = Diagonal(eta .+ D) - D1betaD1 + conds[i] = cond(Matrix(hmD1betaD1)) + end + plot!(Ns, conds, label = "p = $accuracy_order", linewidth = 2, linestyle = :auto) + end + savefig(joinpath(OUT_LAKEATREST, "condition_numbers.pdf")) +end + +# Chapter 4.3 Dingemans +const OUT_DINGEMANS = joinpath(OUT, "dingemans") +ispath(OUT_DINGEMANS) || mkpath(OUT_DINGEMANS) + +# Plot of total waterheight for different models at different points in time +function fig_9() + linewidth = 3 + fontsize = 16 + linestyles = [:solid, :dash, :dot] + + N = 512 + steps = [100, 200, 300, 500] + xlims_zoom = [(-25, 0), (5, 30), (20, 45), (-100, -75)] + ylim_zoom = (0.75, 0.85) + + trixi_include(joinpath(EXAMPLES_DIR_BBMBBM_VARIABLE, + "bbm_bbm_variable_bathymetry_1d_dingemans.jl"); + N = N) + plot(layout = (2, 2), ylim = (-0.05, 0.86), size = (1200, 800), + titlefontsize = fontsize) + for (i, step) in enumerate(steps) + plot!(semi => sol, step = step, conversion = waterheight_total, label = "BBM-BBM", + subplot = i, plot_title = "", linewidth = linewidth, legend = :none, + guidefontsize = fontsize, tickfontsize = fontsize, linestyle = linestyles[1]) + plot!(semi => sol, step = step, inset = (i, bbox(0.1, 0.2, 0.6, 0.5)), + conversion = waterheight_total, linewidth = linewidth, legend = :none, + framestyle = :box, xlim = xlims_zoom[i], ylim = ylim_zoom, + subplot = length(steps) + i, plot_title = "", title = "", xguide = "", + yguide = "", linestyle = linestyles[1]) + end + + trixi_include(joinpath(EXAMPLES_DIR_SVAERD_KALISCH, "svaerd_kalisch_1d_dingemans.jl"); + N = N) + for (i, step) in enumerate(steps) + plot!(semi => sol, step = step, conversion = waterheight_total, + label = "Svärd-Kalisch (set 3)", subplot = i, plot_bathymetry = false, + plot_title = "", linewidth = linewidth, legend = :none, + guidefontsize = fontsize, tickfontsize = fontsize, color = 2, + linestyle = linestyles[2]) + plot!(semi => sol, step = step, conversion = waterheight_total, + linewidth = linewidth, legend = :none, framestyle = :box, + xlim = xlims_zoom[i], ylim = ylim_zoom, subplot = length(steps) + i, + plot_title = "", title = "", xguide = "", yguide = "", color = 2, + linestyle = linestyles[2]) + end + + include("elixir_shallowwater_1d_dingemans.jl") + for (i, step) in enumerate(steps) + plot!(PlotData1D(sol.u[step], semi)["H"], label = "Shallow water", subplot = i, + title = "t = $(round(sol.t[step], digits = 2))", plot_title = "", + linewidth = linewidth, legend = :none, guidefontsize = fontsize, + tickfontsize = fontsize, color = 3, linestyle = linestyles[3]) + plot!(PlotData1D(sol.u[step], semi)["H"], linewidth = linewidth, legend = :none, + framestyle = :box, xlim = xlims_zoom[i], ylim = ylim_zoom, + subplot = length(steps) + i, plot_title = "", title = "", xguide = "", + yguide = "", color = 3, linestyle = linestyles[3]) + end + + # dirty hack to have one legend for all subplots + plot!(subplot = 3, legend_column = 2, bottom_margin = 22 * Plots.mm, + legend = (0.7, -0.34), legendfontsize = 12) + plot!(left_margin = 5 * Plots.mm) + + # plot boxes + for i in 1:length(steps) + plot!([xlims_zoom[i][1], xlims_zoom[i][2]], [ylim_zoom[1], ylim_zoom[1]], + color = :black, label = :none, subplot = i, linewidth = 2) + plot!([xlims_zoom[i][1], xlims_zoom[i][2]], [ylim_zoom[2], ylim_zoom[2]], + color = :black, label = :none, subplot = i, linewidth = 2) + plot!([xlims_zoom[i][1], xlims_zoom[i][1]], [ylim_zoom[1], ylim_zoom[2]], + color = :black, label = :none, subplot = i, linewidth = 2) + plot!([xlims_zoom[i][2], xlims_zoom[i][2]], [ylim_zoom[1], ylim_zoom[2]], + color = :black, label = :none, subplot = i, linewidth = 2) + end + # plot connecting lines + upper_corners = [[-119.5, 0.68], [-9.5, 0.68]] + for i in 1:length(steps) + plot!([xlims_zoom[i][1], upper_corners[1][1]], [ylim_zoom[1], upper_corners[1][2]], + color = :black, label = :none, subplot = i, linewidth = 2) + plot!([xlims_zoom[i][2], upper_corners[2][1]], [ylim_zoom[1], upper_corners[2][2]], + color = :black, label = :none, subplot = i, linewidth = 2) + end + savefig(joinpath(OUT_DINGEMANS, "waterheight_over_time.pdf")) +end + +# Plot of total waterheight for Svärd-Kalisch equations at different points in space and different orders of accuracy +function fig_10() + ylim = (0.75, 0.85) + yticks = [0.76, 0.78, 0.8, 0.82, 0.84] + x_values = [3.04, 9.44, 20.04, 26.04, 30.44, 37.04] + tlims = [ + (15.0, 45.0), + (19.0, 48.0), + (25.0, 52.0), + (30.0, 60.0), + (33.0, 61.0), + (35.0, 65.0), + ] + plot(layout = (3, 2)) + + N = 512 + tspan = (0.0, 70.0) + saveat = range(tspan..., length = 1000) + accuracy_orders = [2, 4, 6] + linestyles = [:solid, :dash, :dot] + + for (i, accuracy_order) in enumerate(accuracy_orders) + trixi_include(joinpath(EXAMPLES_DIR_SVAERD_KALISCH, + "svaerd_kalisch_1d_dingemans.jl"); + N = N, tspan = tspan, accuracy_order = accuracy_order, + saveat = saveat) + for (j, x) in enumerate(x_values) + index = argmin(abs.(DispersiveShallowWater.grid(semi) .- x)) + title = "x = $(round(DispersiveShallowWater.grid(semi)[index], digits = 4))" + plot!(semi => sol, x, conversion = waterheight_total, subplot = j, + xlim = tlims[j], ylim = ylim, plot_title = "", title = title, + legend = nothing, yticks = yticks, linewidth = 2, titlefontsize = 10, + label = "p = $accuracy_order ", linestyle = linestyles[i]) + end + end + + plot!(subplot = 5, legend = (0.82, -1.0), legend_column = 3, legendfontsize = 8, + bottom_margin = 10 * Plots.mm) + savefig(joinpath(OUT_DINGEMANS, "waterheight_at_x_accuracy_order.pdf")) +end + +# Plots of total waterheight for Svärd-Kalisch equations at different points in space and different types of solvers +function fig_11() + ylim = (0.75, 0.85) + yticks = [0.76, 0.78, 0.8, 0.82, 0.84] + x_values = [3.04, 9.44, 20.04, 26.04, 30.44, 37.04] + tlims = [ + (15.0, 45.0), + (19.0, 48.0), + (25.0, 52.0), + (30.0, 60.0), + (33.0, 61.0), + (35.0, 65.0), + ] + plot(layout = (3, 2)) + + N = 512 + tspan = (0.0, 70.0) + saveat = range(tspan..., length = 1000) + accuracy_order = 4 + linestyles = [:solid, :dash, :dot] + + coordinates_min = -138.0 + coordinates_max = 46.0 + p = 3 # N needs to be divisible by p + 1 + D_legendre = legendre_derivative_operator(-1.0, 1.0, p + 1) + uniform_mesh = UniformPeriodicMesh1D(coordinates_min, coordinates_max, div(N, p + 1)) + D1 = couple_discontinuously(D_legendre, uniform_mesh) + D_pl = couple_discontinuously(D_legendre, uniform_mesh, Val(:plus)) + D_min = couple_discontinuously(D_legendre, uniform_mesh, Val(:minus)) + D2 = sparse(D_pl) * sparse(D_min) + solver_DG = Solver(D1, D2) + + p = 4 # N needs to be divisible by p + D_legendre = legendre_derivative_operator(-1.0, 1.0, p + 1) + uniform_mesh = UniformPeriodicMesh1D(coordinates_min, coordinates_max, div(N, p)) + D1 = couple_continuously(D_legendre, uniform_mesh) + D2_legendre = legendre_second_derivative_operator(-1.0, 1.0, p + 1) + D2 = couple_continuously(D2_legendre, uniform_mesh) + solver_CG = Solver(D1, D2) + + solvers = [solver_DG, :none, solver_CG] + labels = ["DG ", "FD ", "CG "] + + for (i, solver) in enumerate(solvers) + if solver == :none + trixi_include(joinpath(EXAMPLES_DIR_SVAERD_KALISCH, + "svaerd_kalisch_1d_dingemans.jl"); + N = N, tspan = tspan, accuracy_order = accuracy_order, + saveat = saveat) + else + trixi_include(joinpath(EXAMPLES_DIR_SVAERD_KALISCH, + "svaerd_kalisch_1d_dingemans.jl"); + N = N, tspan = tspan, accuracy_order = accuracy_order, + saveat = saveat, solver = solvers[i]) + end + for (j, x) in enumerate(x_values) + index = argmin(abs.(DispersiveShallowWater.grid(semi) .- x)) + title = "x = $(round(DispersiveShallowWater.grid(semi)[index], digits = 4))" + plot!(semi => sol, x, conversion = waterheight_total, subplot = j, + xlim = tlims[j], ylim = ylim, plot_title = "", title = title, + legend = nothing, yticks = yticks, linewidth = 2, titlefontsize = 10, + label = labels[i], linestyle = linestyles[i]) + end + end + + plot!(subplot = 5, legend = (0.86, -1.0), legend_column = 3, legendfontsize = 8, + bottom_margin = 10 * Plots.mm) + savefig(joinpath(OUT_DINGEMANS, "waterheight_at_x_solver_types.pdf")) +end + +# Plot solution at different points in space and invariants for entropy conservative and dissipative schemes +function fig_12_13() + ylim = (0.75, 0.85) + yticks = [0.76, 0.78, 0.8, 0.82, 0.84] + x_values = [3.04, 9.44, 20.04, 26.04, 30.44, 37.04] + tlims = [ + (15.0, 45.0), + (19.0, 48.0), + (25.0, 52.0), + (30.0, 60.0), + (33.0, 61.0), + (35.0, 65.0), + ] + p1 = plot(layout = (3, 2)) + + N = 512 + tspan = (0.0, 70.0) + saveat = range(tspan..., length = 1000) + accuracy_order = 4 + linestyles = [:solid, :dash, :dot] + linewidth = 2 + titlefontsize = 10 + + labels = ["EC baseline ", "EC relaxation ", "ED "] + + function plot_at_x(semi, sol, i) + for (j, x) in enumerate(x_values) + index = argmin(abs.(DispersiveShallowWater.grid(semi) .- x)) + title = "x = $(round(DispersiveShallowWater.grid(semi)[index], digits = 4))" + plot!(p1, semi => sol, x, conversion = waterheight_total, subplot = j, + xlim = tlims[j], ylim = ylim, plot_title = "", title = title, + legend = nothing, yticks = yticks, linewidth = linewidth, + titlefontsize = titlefontsize, + label = labels[i], linestyle = linestyles[i]) + end + end + + trixi_include(joinpath(EXAMPLES_DIR_SVAERD_KALISCH, "svaerd_kalisch_1d_dingemans.jl"); + N = N, tspan = tspan, accuracy_order = accuracy_order, saveat = saveat) + plot_at_x(semi, sol, 1) + p2 = plot(analysis_callback, title = labels[1], legend = :none, + linestyles = [:solid :dash :dot], + linewidth = linewidth, layout = 3, subplot = 1, titlefontsize = titlefontsize) + + trixi_include(joinpath(EXAMPLES_DIR_SVAERD_KALISCH, + "svaerd_kalisch_1d_dingemans_relaxation.jl"); + N = N, tspan = tspan, accuracy_order = accuracy_order, saveat = saveat) + plot_at_x(semi, sol, 2) + plot!(p2, analysis_callback, title = labels[2], legend = :none, + linestyles = [:solid :dash :dot], + linewidth = linewidth, subplot = 2, titlefontsize = titlefontsize) + + trixi_include(joinpath(EXAMPLES_DIR_SVAERD_KALISCH, + "svaerd_kalisch_1d_dingemans_upwind.jl"); + N = N, tspan = tspan, accuracy_order = accuracy_order, saveat = saveat) + plot_at_x(semi, sol, 3) + plot!(p2, analysis_callback, title = labels[3], legend = :none, + linestyles = [:solid :dash :dot], + linewidth = linewidth, subplot = 3, titlefontsize = titlefontsize) + + plot!(p1, subplot = 5, legend = (0.86, -1.0), legend_column = 3, legendfontsize = 8, + bottom_margin = 10 * Plots.mm) + plot!(p2, subplot = 3, legend = (1.3, 0.6), legendfontsize = 8) + savefig(p1, joinpath(OUT_DINGEMANS, "waterheight_at_x_ec.pdf")) + savefig(p2, joinpath(OUT_DINGEMANS, "invariants_ec.pdf")) +end + +fig_1() +fig_2() +fig_3() +fig_4_5_6() +fig_7() +# fig_orders_different_stencils() +fig_8() +# fig_condition_number() +fig_9() +fig_10() +fig_11() +fig_12_13() diff --git a/visualization/elixir_shallowwater_1d_dingemans.jl b/visualization/elixir_shallowwater_1d_dingemans.jl new file mode 100644 index 00000000..87e83522 --- /dev/null +++ b/visualization/elixir_shallowwater_1d_dingemans.jl @@ -0,0 +1,59 @@ +using Trixi +using OrdinaryDiffEq + +equations = ShallowWaterEquations1D(gravity_constant = 9.81, H0 = 0.8) + +function initial_condition_dingemans_trixi(x, t, equations::ShallowWaterEquations1D) + eta0 = 0.8 + A = 0.02 + # omega = 2*pi/(2.02*sqrt(2)) + k = 0.8406220896381442 # precomputed result of find_zero(k -> omega^2 - equations.gravity * k * tanh(k * eta0), 1.0) using Roots.jl + if x[1] < -30.5 * pi / k || x[1] > -8.5 * pi / k + h = 0.0 + else + h = A * cos(k * x[1]) + end + v = sqrt(equations.gravity / k * tanh(k * eta0)) * h / eta0 + if x[1] < 11.01 || x[1] >= 33.07 + b = 0.0 + elseif 11.01 <= x[1] && x[1] < 23.04 + b = 0.6 * (x[1] - 11.01) / (23.04 - 11.01) + elseif 23.04 <= x[1] && x[1] < 27.04 + b = 0.6 + elseif 27.04 <= x[1] && x[1] < 33.07 + b = 0.6 * (33.07 - x[1]) / (33.07 - 27.04) + else + error("should not happen") + end + eta = h + eta0 + D = -b + return Trixi.prim2cons(SVector(eta, v, b), equations) +end + +initial_condition = initial_condition_dingemans_trixi + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal) +solver = DGSEM(polydeg = 3, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = -138.0 +coordinates_max = 46.0 +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 7, + n_cells_max = 10_000) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) +tspan = (0.0, 70.0) +ode = Trixi.semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_callback = Trixi.AnalysisCallback(semi, interval = 100) + +callbacks = CallbackSet(summary_callback, analysis_callback) + +saveat = range(tspan..., length = 500) +sol = solve(ode, Tsit5(), reltol = 1e-7, abstol = 1e-7, + save_everystep = false, callback = callbacks, saveat = saveat); +summary_callback() # print the timer summary diff --git a/plot_examples.jl b/visualization/plot_examples.jl similarity index 99% rename from plot_examples.jl rename to visualization/plot_examples.jl index 3fb4e6fe..58e2adb8 100644 --- a/plot_examples.jl +++ b/visualization/plot_examples.jl @@ -1,4 +1,4 @@ -# This script solves all the examples in `examples_dir()` and creates a .gif of solution +# This script solves all the examples in `examples_dir()` and creates a .gif of the solution # as well as some additional plots. All plots are saved in a directory `out/` and # subdirectories therein.