diff --git a/.github/workflows/Format-check.yml b/.github/workflows/Format-check.yml index 41987435..9edf4ae6 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"], verbose = true)' + julia -e 'using JuliaFormatter; format(["src", "test", "examples", "create_figures.jl", "plot_examples.jl"], verbose = true)' - name: Format check run: | julia -e ' diff --git a/create_figures.jl b/create_figures.jl index c4d2a2a5..272fa3c1 100644 --- a/create_figures.jl +++ b/create_figures.jl @@ -17,17 +17,23 @@ function fig_1() 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) + 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) + 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) + 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) + 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 @@ -41,28 +47,32 @@ function fig_2() g = 1.0 c0 = sqrt(g * h0) - k = 0.01:0.5:(8*pi) + 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, ylims = (0.54, 1.0), inset = bbox(0.35, 0.1, 0.35, 0.3), subplot = 2, legend = nothing, linewidth = linewidth, markershape = :circle, markersize = markersize) + 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), + inset = bbox(0.35, 0.1, 0.35, 0.3), subplot = 2, legend = nothing, + linewidth = linewidth, markershape = :circle, markersize = markersize) 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) + 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, 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 + alpha_set1 = -1 / 3 * c0 * h0^2 beta_set1 = 0.0 * h0^3 gamma_set1 = 0.0 * c0 * h0^3 @@ -79,11 +89,11 @@ function fig_2() 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) + 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 @@ -112,24 +122,30 @@ function fig_3() 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") + plot([], label = :none, xscale = :log2, yscale = :log10, xticks = all_Ns, xlabel = "N", + ylabel = L"||\eta - \eta_{ana}||_2 + ||v - v_{ana}||_2") 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]) + _, 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) + 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") end savefig(joinpath(OUT, "orders.pdf")) - plot([], label = :none, xscale = :log2, yscale = :log10, xticks = all_Ns, xlabel = "N", ylabel = L"||\eta - \eta_{ana}||_2 + ||v - v_{ana}||_2") + plot([], label = :none, xscale = :log2, yscale = :log10, xticks = all_Ns, xlabel = "N", + ylabel = L"||\eta - \eta_{ana}||_2 + ||v - v_{ana}||_2") 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]) + _, 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) + 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") diff --git a/plot_examples.jl b/plot_examples.jl index dd03b6c3..00c8e606 100644 --- a/plot_examples.jl +++ b/plot_examples.jl @@ -66,38 +66,38 @@ 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)], - tspan=(0.0, 30.0)) + ylims_gif=[(-8, 4), (-10, 30)], + tspan=(0.0, 30.0)) ############################################################################### # Travelling wave solution for one-dimensional BBM-BBM equations with periodic boundary conditions # using periodic SBP operators. Uses the BBM-BBM equations with variable bathymetry, but sets the bathymetry # 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], - tspan=(0.0, 50.0)) + "bbm_bbm_variable_bathymetry_1d_basic.jl"), + ylims_gif=[(-8, 4), :auto], + tspan=(0.0, 50.0)) ############################################################################### # One-dimensional BBM-BBM equations with a Gaussian bump as initial condition for the water height # and initially still water. The bathymetry is a sine function. Relaxation is used, so the solution # 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)], - tspan=(0.0, 10.0)) + "bbm_bbm_variable_bathymetry_1d_relaxation.jl"), + ylims_gif=[(-1.5, 6.0), (-10.0, 10.0)], + tspan=(0.0, 10.0)) ############################################################################### # One-dimensional BBM-BBM equations with a Gaussian bump as initial condition for the water height @@ -113,9 +113,9 @@ const EXAMPLES_DIR_SVAERD_KALISCH = "svaerd_kalisch_1d" # and initially still water. The bathymetry is a sine function. Relaxation is used, so the solution # 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)], - tspan=(0.0, 10.0)) + "bbm_bbm_variable_bathymetry_1d_upwind_relaxation.jl"), + ylims_gif=[(-1.5, 6.0), (-10.0, 10.0)], + tspan=(0.0, 10.0)) ############################################################################### # One-dimensional BBM-BBM equations with a constant water height @@ -123,82 +123,82 @@ const EXAMPLES_DIR_SVAERD_KALISCH = "svaerd_kalisch_1d" # is energy-conservative. Uses periodic finite difference SBP operators. The solution should be # (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)], - tspan=(0.0, 10.0)) + "bbm_bbm_variable_bathymetry_1d_well_balanced.jl"), + ylims_gif=[(2.0 - 1e-3, 2.0 + 1e-3), (-1e-3, 1e-3)], + tspan=(0.0, 10.0)) ############################################################################### # One-dimensional BBM-BBM equations with initial condition that models # a wave make. This setup comes from experiments by W. M. Dingemans. @plot_example(joinpath(EXAMPLES_DIR_BBMBBM_VARIABLE, - "bbm_bbm_variable_bathymetry_1d_dingemans.jl"), - ylims_gif=[(-0.1, 0.9), (-0.3, 0.3)], - ylims_x=[:auto, :auto], - 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), - ], - tspan=(0.0, 70.0)) + "bbm_bbm_variable_bathymetry_1d_dingemans.jl"), + 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), + (19.0, 48.0), + (25.0, 52.0), + (30.0, 60.0), + (33.0, 61.0), + (35.0, 65.0), + ], + tspan=(0.0, 70.0)) ############################################################################### # One-dimensional equations from Svärd and Kalisch with initial condition that models # 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], - 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), - ], - tspan=(0.0, 70.0)) + "svaerd_kalisch_1d_dingemans.jl"), + 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), + (19.0, 48.0), + (25.0, 52.0), + (30.0, 60.0), + (33.0, 61.0), + (35.0, 65.0), + ], + tspan=(0.0, 70.0)) ############################################################################### # One-dimensional equations from Svärd and Kalisch with initial condition that models # 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], - 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), - ], - tspan=(0.0, 70.0)) + "svaerd_kalisch_1d_dingemans_upwind.jl"), + 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), + (19.0, 48.0), + (25.0, 52.0), + (30.0, 60.0), + (33.0, 61.0), + (35.0, 65.0), + ], + tspan=(0.0, 70.0)) ############################################################################### # One-dimensional equations from Svärd and Kalisch with initial condition that models # a wave make. This setup comes from experiments by W. M. Dingemans. Relaxation is used # 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], - 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), - ], - tspan=(0.0, 70.0)) + "svaerd_kalisch_1d_dingemans_relaxation.jl"), + 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), + (19.0, 48.0), + (25.0, 52.0), + (30.0, 60.0), + (33.0, 61.0), + (35.0, 65.0), + ], + tspan=(0.0, 70.0)) ############################################################################### # One-dimensional Svärd-Kalisch equations with a constant water height @@ -206,6 +206,6 @@ const EXAMPLES_DIR_SVAERD_KALISCH = "svaerd_kalisch_1d" # is energy-conservative. Uses periodic finite difference SBP operators. The solution should be # (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)], - tspan=(0.0, 10.0)) + "svaerd_kalisch_1d_well_balanced.jl"), + ylims=[(2.0 - 1e-3, 2.0 + 1e-3), (-1e-3, 1e-3)], + tspan=(0.0, 10.0)) diff --git a/src/util.jl b/src/util.jl index 5cc41d55..c4267ac5 100644 --- a/src/util.jl +++ b/src/util.jl @@ -193,7 +193,7 @@ function analyze_convergence(errors, iterations, # For the following iterations print errors and EOCs for j in 2:iterations for k in 1:nvariables - @printf("%-5d", initial_N * 2^(j - 1)) + @printf("%-5d", initial_N*2^(j - 1)) @printf("%-10.2e", error[j, k]) @printf("%-10.2f", eocs[kind][j - 1, k]) end