From a1a0fb38035657b225271cdec0404813fdb8350e Mon Sep 17 00:00:00 2001 From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Date: Mon, 25 Sep 2023 01:16:30 +0200 Subject: [PATCH] Add `convergence_test` (#43) * add convergece_test * format * explicitly define accuracy_order in convergence_test relaxation example --- create_figures.jl | 176 +++++++++---------- examples/bbm_bbm_1d/bbm_bbm_1d_basic.jl | 3 +- examples/bbm_bbm_1d/bbm_bbm_1d_relaxation.jl | 3 +- src/DispersiveShallowWater.jl | 4 +- src/callbacks_step/analysis.jl | 12 ++ src/util.jl | 135 +++++++++++++- test/test_util.jl | 14 +- 7 files changed, 245 insertions(+), 102 deletions(-) diff --git a/create_figures.jl b/create_figures.jl index 85440529..ff4f23ce 100644 --- a/create_figures.jl +++ b/create_figures.jl @@ -2,8 +2,8 @@ using DispersiveShallowWater using Plots # Use a macro to avoid world age issues when defining new initial conditions etc. -# inside an elixir. -macro plot_elixir(filename, args...) +# inside an example. +macro plot_example(filename, args...) local ylims_gif = get_kwarg(args, :ylims_gif, nothing) local ylims_x = get_kwarg(args, :ylims_x, nothing) local x_values = get_kwarg(args, :x_values, []) @@ -54,46 +54,46 @@ function get_kwarg(args, keyword, default_value) return val end -EXAMPLES_DIR_BBMBBM = "bbm_bbm_1d" -EXAMPLES_DIR_BBMBBM_VARIABLE = "bbm_bbm_variable_bathymetry_1d" -EXAMPLES_DIR_SVAERD_KALISCH = "svaerd_kalisch_1d" +const EXAMPLES_DIR_BBMBBM = "bbm_bbm_1d" +const EXAMPLES_DIR_BBMBBM_VARIABLE = "bbm_bbm_variable_bathymetry_1d" +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_elixir(joinpath(EXAMPLES_DIR_BBMBBM, "bbm_bbm_1d_basic.jl"), - ylims_gif=[(-8, 4), :auto], tspan=(0.0, 50.0)) +@plot_example(joinpath(EXAMPLES_DIR_BBMBBM, "bbm_bbm_1d_basic.jl"), + 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_elixir(joinpath(EXAMPLES_DIR_BBMBBM, "bbm_bbm_1d_dg.jl"), - ylims_gif=[(-4, 2), :auto]) +@plot_example(joinpath(EXAMPLES_DIR_BBMBBM, "bbm_bbm_1d_dg.jl"), + 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_elixir(joinpath(EXAMPLES_DIR_BBMBBM, "bbm_bbm_1d_relaxation.jl"), - ylims_gif=[(-8, 4), (-10, 30)], - tspan=(0.0, 30.0)) +@plot_example(joinpath(EXAMPLES_DIR_BBMBBM, "bbm_bbm_1d_relaxation.jl"), + 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_elixir(joinpath(EXAMPLES_DIR_BBMBBM_VARIABLE, - "bbm_bbm_variable_bathymetry_1d_basic.jl"), - ylims_gif=[(-8, 4), :auto], - tspan=(0.0, 50.0)) +@plot_example(joinpath(EXAMPLES_DIR_BBMBBM_VARIABLE, + "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_elixir(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)) +@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)) ############################################################################### # One-dimensional BBM-BBM equations with a Gaussian bump as initial condition for the water height @@ -108,100 +108,100 @@ EXAMPLES_DIR_SVAERD_KALISCH = "svaerd_kalisch_1d" # 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 discontinuously coupled SBP operators. -@plot_elixir(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)) +@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)) ############################################################################### # One-dimensional BBM-BBM equations with a constant water height # and initially still water. The bathymetry is discontinuous. Relaxation is used, so the solution # is energy-conservative. Uses periodic finite difference SBP operators. The solution should be # (exactly) constant in time. -@plot_elixir(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)) +@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)) ############################################################################### # One-dimensional BBM-BBM equations with initial condition that models # a wave make. This setup comes from experiments by W. M. Dingemans. -@plot_elixir(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)) +@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)) ############################################################################### # 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_elixir(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)) +@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)) ############################################################################### # 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_elixir(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)) +@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)) ############################################################################### # 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_elixir(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)) +@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)) ############################################################################### # One-dimensional Svärd-Kalisch equations with a constant water height # and initially still water. The bathymetry is discontinuous. Relaxation is used, so the solution # is energy-conservative. Uses periodic finite difference SBP operators. The solution should be # (exactly) constant in time. -@plot_elixir(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)) +@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)) diff --git a/examples/bbm_bbm_1d/bbm_bbm_1d_basic.jl b/examples/bbm_bbm_1d/bbm_bbm_1d_basic.jl index 234ac9b9..00ac9e29 100644 --- a/examples/bbm_bbm_1d/bbm_bbm_1d_basic.jl +++ b/examples/bbm_bbm_1d/bbm_bbm_1d_basic.jl @@ -17,7 +17,8 @@ N = 512 mesh = Mesh1D(coordinates_min, coordinates_max, N + 1) # 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_1d/bbm_bbm_1d_relaxation.jl b/examples/bbm_bbm_1d/bbm_bbm_1d_relaxation.jl index 3062d37d..b5c49100 100644 --- a/examples/bbm_bbm_1d/bbm_bbm_1d_relaxation.jl +++ b/examples/bbm_bbm_1d/bbm_bbm_1d_relaxation.jl @@ -17,7 +17,8 @@ N = 512 mesh = Mesh1D(coordinates_min, coordinates_max, N + 1) # 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 620b63bd..5757bb2b 100644 --- a/src/DispersiveShallowWater.jl +++ b/src/DispersiveShallowWater.jl @@ -20,7 +20,6 @@ using SummationByPartsOperators: AbstractDerivativeOperator, derivative_order, integrate import SummationByPartsOperators: grid, xmin, xmax -include("util.jl") include("boundary_conditions.jl") include("mesh.jl") include("solver.jl") @@ -28,8 +27,9 @@ include("equations/equations.jl") include("semidiscretization.jl") include("callbacks_step/callbacks_step.jl") include("visualization.jl") +include("util.jl") -export examples_dir, get_examples, default_example, trixi_include +export examples_dir, get_examples, default_example, trixi_include, convergence_test export BBMBBMEquations1D, BBMBBMVariableEquations1D, SvärdKalischEquations1D, SvaerdKalischEquations1D diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index 93806d37..8ff8e6da 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -307,6 +307,18 @@ function analyze_integrals!(current_integrals, i, analysis_integrals::Tuple{}, u nothing end +# used for error checks and EOC analysis +function (cb::DiscreteCallback{Condition, Affect!})(sol) where {Condition, + Affect! <: + AnalysisCallback} + analysis_callback = cb.affect! + semi = sol.prob.p + + l2_error, linf_error = calc_error_norms(sol.u[end], sol.t[end], semi) + + return (; l2 = l2_error, linf = linf_error) +end + function analyze(quantity, u_ode, t, semi::Semidiscretization) integrate_quantity(u_ode -> quantity(u_ode, semi.equations), u_ode, semi) end diff --git a/src/util.jl b/src/util.jl index fa9d0aa2..9e822f21 100644 --- a/src/util.jl +++ b/src/util.jl @@ -95,11 +95,127 @@ julia> redirect_stdout(devnull) do 0.1 ``` """ -function trixi_include(mod::Module, elixir::AbstractString; kwargs...) - Base.include(ex -> replace_assignments(insert_maxiters(ex); kwargs...), mod, elixir) +function trixi_include(mod::Module, example::AbstractString; kwargs...) + Base.include(ex -> replace_assignments(insert_maxiters(ex); kwargs...), mod, example) end -trixi_include(elixir::AbstractString; kwargs...) = trixi_include(Main, elixir; kwargs...) +trixi_include(example::AbstractString; kwargs...) = trixi_include(Main, example; kwargs...) + +""" + convergence_test([mod::Module=Main,] example::AbstractString, iterations; kwargs...) + +Run `iterations` simulations using the setup given in `example` and compute +the experimental order of convergence (EOC) in the ``L^2`` and ``L^\\infty`` norm. +In each iteration, the resolution of the respective mesh will be doubled. +Additional keyword arguments `kwargs...` and the optional module `mod` are passed directly +to [`trixi_include`](@ref). + +Adjusted from [Trixi.jl](https://github.com/trixi-framework/Trixi.jl). +""" +function convergence_test(mod::Module, example::AbstractString, iterations; kwargs...) + @assert(iterations>1, + "Number of iterations must be bigger than 1 for a convergence analysis") + + # Types of errors to be calculated + errors = Dict(:l2 => Float64[], :linf => Float64[]) + + initial_N = extract_initial_N(example, kwargs) + + # run simulations and extract errors + for iter in 1:iterations + println("Running convtest iteration ", iter, "/", iterations) + + trixi_include(mod, example; kwargs..., N = initial_N * 2^(iter - 1)) + + l2_error, linf_error = mod.analysis_callback(mod.sol) + + # collect errors as one vector to reshape later + append!(errors[:l2], l2_error) + append!(errors[:linf], linf_error) + + println("\n\n") + println("#"^100) + end + + # Use raw error values to compute EOC + analyze_convergence(errors, iterations, mod.semi) +end + +# Analyze convergence for any semidiscretization +# Note: this intermediate method is to allow dispatching on the semidiscretization +function analyze_convergence(errors, iterations, semi::Semidiscretization) + _, equations, _, _ = mesh_equations_solver_cache(semi) + variablenames = varnames(prim2prim, equations) + analyze_convergence(errors, iterations, variablenames) +end + +# This method is called with the collected error values to actually compute and print the EOC +function analyze_convergence(errors, iterations, + variablenames::Union{Tuple, AbstractArray}) + nvariables = length(variablenames) + + # Reshape errors to get a matrix where the i-th row represents the i-th iteration + # and the j-th column represents the j-th variable + errorsmatrix = Dict(kind => transpose(reshape(error, (nvariables, iterations))) + for (kind, error) in errors) + + # Calculate EOCs where the columns represent the variables + # As dx halves in every iteration the denominator needs to be log(1/2) + eocs = Dict(kind => log.(error[2:end, :] ./ error[1:(end - 1), :]) ./ log(1 / 2) + for (kind, error) in errorsmatrix) + + eoc_mean_values = Dict{Symbol, Any}() + eoc_mean_values[:variables] = variablenames + + for (kind, error) in errorsmatrix + println(kind) + + for v in variablenames + @printf("%-20s", v) + end + println("") + + for k in 1:nvariables + @printf("%-10s", "error") + @printf("%-10s", "EOC") + end + println("") + + # Print errors for the first iteration + for k in 1:nvariables + @printf("%-10.2e", error[1, k]) + @printf("%-10s", "-") + end + println("") + + # For the following iterations print errors and EOCs + for j in 2:iterations + for k in 1:nvariables + @printf("%-10.2e", error[j, k]) + @printf("%-10.2f", eocs[kind][j - 1, k]) + end + println("") + end + println("") + + # Print mean EOCs + mean_values = zeros(nvariables) + for v in 1:nvariables + mean_values[v] = sum(eocs[kind][:, v]) ./ length(eocs[kind][:, v]) + @printf("%-10s", "mean") + @printf("%-10.2f", mean_values[v]) + end + eoc_mean_values[kind] = mean_values + println("") + println("-"^100) + end + + return eoc_mean_values, errorsmatrix +end + +function convergence_test(example::AbstractString, iterations; kwargs...) + convergence_test(Main, example::AbstractString, iterations; kwargs...) +end # Helper methods used in the functions defined above, also copied from Trixi.jl @@ -179,6 +295,19 @@ function find_assignment(expr, destination) result end +function extract_initial_N(example, kwargs) + code = read(example, String) + expr = Meta.parse("begin \n$code \nend") + + if haskey(kwargs, :N) + return kwargs[:N] + else + # get N from the example + N = find_assignment(expr, :N) + return N + end +end + """ @autoinfiltrate @autoinfiltrate condition::Bool diff --git a/test/test_util.jl b/test/test_util.jl index 534c4803..8a356847 100644 --- a/test/test_util.jl +++ b/test/test_util.jl @@ -1,9 +1,9 @@ using Test: @test # Use a macro to avoid world age issues when defining new initial conditions etc. -# inside an elixir. +# inside an example. """ - @test_trixi_include(elixir; l2=nothing, linf=nothing, cons_error=nothing + @test_trixi_include(example; l2=nothing, linf=nothing, cons_error=nothing change_waterheight=nothing, change_velocity=nothing, change_entropy=nothing, @@ -12,7 +12,7 @@ using Test: @test atol=1e-12, rtol=sqrt(eps()), atol_ints=1e-11, rtol_ints=sqrt(eps())) -Test by calling `trixi_include(elixir; parameters...)`. +Test by calling `trixi_include(example; parameters...)`. By default, only the absence of error output is checked. If `l2`, `linf` or `cons_error` are specified, in addition the resulting L2/Linf/conservation errors are compared approximately against these reference values, using `atol, rtol` @@ -21,7 +21,7 @@ If `change_waterheight`, `change_velocity`, `change_momemtum`, `change_entropy`, or `lake_at_rest` are specified, in addition the resulting changes of the different errors are compared approximately against these reference values, using `atol_ints`, `rtol_ints` as absolute/relative tolerance. """ -macro test_trixi_include(elixir, args...) +macro test_trixi_include(example, args...) local l2 = get_kwarg(args, :l2, nothing) local linf = get_kwarg(args, :linf, nothing) local cons_error = get_kwarg(args, :cons_error, nothing) @@ -49,10 +49,10 @@ macro test_trixi_include(elixir, args...) quote println("═"^100) - println($elixir) + println($example) # evaluate examples in the scope of the module they're called from - @test_nowarn trixi_include(@__MODULE__, $elixir; $kwargs...) + @test_nowarn trixi_include(@__MODULE__, $example; $kwargs...) # if present, compare l2, linf and conservation errors against reference values if !isnothing($l2) || !isnothing($linf) || !isnothing($cons_error) @@ -173,7 +173,7 @@ macro trixi_testset(name, expr) using DispersiveShallowWater include(@__FILE__) # We define `EXAMPLES_DIR` in (nearly) all test modules and use it to - # get the path to the elixirs to be tested. However, that's not required + # get the path to the examples to be tested. However, that's not required # and we want to fail gracefully if it's not defined. try import ..EXAMPLES_DIR