diff --git a/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans_cg.jl b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans_cg.jl new file mode 100644 index 00000000..25ffc1a8 --- /dev/null +++ b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans_cg.jl @@ -0,0 +1,47 @@ +using OrdinaryDiffEq +using DispersiveShallowWater +using SummationByPartsOperators: legendre_derivative_operator, UniformPeriodicMesh1D, + couple_continuously, legendre_second_derivative_operator + +############################################################################### +# Semidiscretization of the Svärd-Kalisch equations + +equations = SvaerdKalischEquations1D(gravity_constant = 9.81, eta0 = 0.8, alpha = 0.0, + beta = 0.27946992481203003, gamma = 0.0521077694235589) + +initial_condition = initial_condition_dingemans +boundary_conditions = boundary_condition_periodic + +# create homogeneous mesh +coordinates_min = -138.0 +coordinates_max = 46.0 +N = 512 +mesh = Mesh1D(coordinates_min, coordinates_max, N) + +# create solver with periodic SBP operators of accuracy order 4 +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 = 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, 70.0) +ode = semidiscretize(semi, tspan) +analysis_callback = AnalysisCallback(semi; interval = 10, + extra_analysis_errors = (:conservation_error,), + extra_analysis_integrals = (waterheight_total, + entropy, + entropy_modified)) +callbacks = CallbackSet(analysis_callback) + +saveat = range(tspan..., length = 500) +sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7, + save_everystep = false, callback = callbacks, saveat = saveat) diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index cc54d153..4c8ac038 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -55,7 +55,7 @@ function Base.show(io::IO, ::MIME"text/plain", println(io, " error ", idx, ": ", error) end for (idx, integral) in enumerate(analysis_integrals) - println(io, " integral ", idx, ": ", integral) + println(io, " integral ", idx, ": ", pretty_form_utf(integral)) end end end diff --git a/test/Project.toml b/test/Project.toml index 66b4d397..d699460a 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,6 +1,7 @@ [deps] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SummationByPartsOperators = "9f78cca6-572e-554e-b819-917d2f1cf240" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" @@ -8,4 +9,5 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] Aqua = "0.7" OrdinaryDiffEq = "6.49.1" +Plots = "1.16" SummationByPartsOperators = "0.5.41" diff --git a/test/runtests.jl b/test/runtests.jl index a3bc94a7..55376424 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,6 +3,7 @@ using Test @testset "DispersiveShallowWater.jl" begin include("test_aqua.jl") include("test_unit.jl") + include("test_visualization.jl") include("test_bbm_bbm_1d.jl") include("test_bbm_bbm_variable_bathymetry_1d.jl") include("test_svaerd_kalisch_1d.jl") diff --git a/test/test_svaerd_kalisch_1d.jl b/test/test_svaerd_kalisch_1d.jl index aec84ee6..935e22b4 100644 --- a/test/test_svaerd_kalisch_1d.jl +++ b/test/test_svaerd_kalisch_1d.jl @@ -21,6 +21,18 @@ EXAMPLES_DIR = joinpath(examples_dir(), "svaerd_kalisch_1d") change_entropy_modified=3.240868750253867e-6) end + @trixi_testset "svaerd_kalisch_1d_dingemans_cg" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "svaerd_kalisch_1d_dingemans_cg.jl"), + tspan=(0.0, 1.0), + l2=[0.22798490823942433 0.7520004851600044 0.0], + linf=[0.03673010870720128 0.12074632168110239 0.0], + cons_error=[1.4210854715202004e-13 4.953054817174909e-5 0.0], + change_waterheight=-1.4210854715202004e-13, + change_entropy=-0.0002425303440531934, + change_entropy_modified=2.9729466177741415e-6) + end + @trixi_testset "svaerd_kalisch_1d_dingemans_upwind" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "svaerd_kalisch_1d_dingemans_upwind.jl"), diff --git a/test/test_unit.jl b/test/test_unit.jl index 80b0acaf..ebe64afc 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -7,31 +7,35 @@ using SummationByPartsOperators: PeriodicDerivativeOperator, UniformPeriodicCoup using SummationByPartsOperators: derivative_order, periodic_derivative_operator, legendre_derivative_operator, UniformPeriodicMesh1D, couple_discontinuously, - upwind_operators + upwind_operators, couple_continuously, + legendre_second_derivative_operator + using SparseArrays: sparse, SparseMatrixCSC @testset "Unit tests" begin @testset "Mesh1D" begin mesh = @test_nowarn Mesh1D(-1, 1, 10) + @test_nowarn print(mesh) + @test_nowarn display(mesh) @test ndims(mesh) == 1 @test xmin(mesh) == -1 @test xmax(mesh) == 1 @test nnodes(mesh) == 10 @test real(mesh) == Int64 - @test_nowarn show(stdout, mesh) end @testset "Solver" begin mesh = Mesh1D(-1.0, 1.0, 10) p = 3 solver = @test_nowarn Solver(mesh, p) + @test_nowarn print(solver) + @test_nowarn display(solver) @test solver.D1 isa PeriodicDerivativeOperator @test solver.D2 isa PeriodicDerivativeOperator @test derivative_order(solver.D1) == 1 @test derivative_order(solver.D2) == 2 @test grid(solver) == grid(solver.D1) == grid(solver.D2) @test real(solver) == Float64 - @test_nowarn show(stdout, solver) D_legendre = legendre_derivative_operator(-1.0, 1.0, p + 1) uniform_mesh = UniformPeriodicMesh1D(-1.0, 1.0, 512 ÷ (p + 1)) @@ -54,10 +58,61 @@ using SparseArrays: sparse, SparseMatrixCSC @test solver.D1 isa PeriodicUpwindOperators @test solver.D2 isa SparseMatrixCSC @test derivative_order(solver.D1) == 1 + + p = 4 # N needs to be divisible by p + D_legendre = legendre_derivative_operator(-1.0, 1.0, p + 1) + uniform_mesh = UniformPeriodicMesh1D(-1.0, 1.0, div(12, 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 = @test_nowarn Solver(D1, D2) + @test solver.D1 isa UniformPeriodicCoupledOperator + @test solver.D2 isa UniformPeriodicCoupledOperator + end + + @testset "Semidiscretization" begin + equations = BBMBBMEquations1D(gravity_constant = 9.81, D = 2.0) + initial_condition = initial_condition_convergence_test + boundary_conditions = boundary_condition_periodic + mesh = Mesh1D(-1, 1, 10) + solver = Solver(mesh, 4) + semi = @test_nowarn Semidiscretization(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + @test_nowarn print(semi) + @test_nowarn display(semi) + @test ndims(semi) == ndims(mesh) == 1 + @test DispersiveShallowWater.eachnode(semi) == DispersiveShallowWater.eachnode(mesh) + @test grid(semi) == grid(solver) + mesh, equations, solver, cache = @test_nowarn DispersiveShallowWater.mesh_equations_solver_cache(semi) + @test mesh == mesh + @test equations == equations + @test solver == solver + end + + @testset "Boundary conditions" begin + boundary_conditions = boundary_condition_periodic + @test_nowarn print(boundary_conditions) + @test_nowarn display(boundary_conditions) end @testset "BBMBBMEquations1D" begin equations = @test_nowarn BBMBBMEquations1D(gravity_constant = 9.81, D = 2.0) + @test_nowarn print(equations) + @test_nowarn display(equations) + conversion_functions = [ + waterheight_total, + waterheight, + velocity, + momentum, + discharge, + entropy, + energy_total, + prim2cons, + prim2prim, + ] + for conversion in conversion_functions + @test DispersiveShallowWater.varnames(conversion, equations) isa Tuple + end q = [42.0, 2.0] @test prim2prim(q, equations) == q @test isapprox(cons2prim(prim2cons(q, equations), equations), q) @@ -67,11 +122,26 @@ using SparseArrays: sparse, SparseMatrixCSC @test momentum(q, equations) == 88.0 @test discharge(q, equations) == 88.0 @test isapprox(energy_total(q, equations), 8740.42) - @test_nowarn show(stdout, equations) end @testset "BBMBBMVariableEquations1D" begin equations = @test_nowarn BBMBBMVariableEquations1D(gravity_constant = 9.81) + @test_nowarn print(equations) + @test_nowarn display(equations) + conversion_functions = [ + waterheight_total, + waterheight, + velocity, + momentum, + discharge, + entropy, + energy_total, + prim2cons, + prim2prim, + ] + for conversion in conversion_functions + @test DispersiveShallowWater.varnames(conversion, equations) isa Tuple + end q = [42.0, 2.0, 2.0] @test prim2prim(q, equations) == q @test isapprox(cons2prim(prim2cons(q, equations), equations), q) @@ -81,7 +151,6 @@ using SparseArrays: sparse, SparseMatrixCSC @test momentum(q, equations) == 88.0 @test discharge(q, equations) == 88.0 @test isapprox(energy_total(q, equations), 8740.42) - @test_nowarn show(stdout, equations) end @testset "SvaerdKalischEquations1D" begin @@ -89,6 +158,24 @@ using SparseArrays: sparse, SparseMatrixCSC alpha = 0.0004040404040404049, beta = 0.49292929292929294, gamma = 0.15707070707070708) + @test_nowarn print(equations) + @test_nowarn display(equations) + conversion_functions = [ + waterheight_total, + waterheight, + velocity, + momentum, + discharge, + entropy, + energy_total, + prim2cons, + prim2prim, + energy_total_modified, + entropy_modified, + ] + for conversion in conversion_functions + @test DispersiveShallowWater.varnames(conversion, equations) isa Tuple + end q = [42.0, 2.0, 2.0] @test prim2prim(q, equations) == q @test isapprox(cons2prim(prim2cons(q, equations), equations), q) @@ -98,7 +185,51 @@ using SparseArrays: sparse, SparseMatrixCSC @test momentum(q, equations) == 88.0 @test discharge(q, equations) == 88.0 @test isapprox(energy_total(q, equations), 8740.42) - @test_nowarn show(stdout, equations) + end + + @testset "AnalysisCallback" begin + equations = SvaerdKalischEquations1D(gravity_constant = 9.81) + initial_condition = initial_condition_dingemans + boundary_conditions = boundary_condition_periodic + mesh = Mesh1D(-1, 1, 10) + solver = Solver(mesh, 4) + semi = Semidiscretization(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + analysis_callback = AnalysisCallback(semi; interval = 10, + extra_analysis_errors = (:conservation_error,), + extra_analysis_integrals = (waterheight_total, + velocity, momentum, + discharge, entropy, + energy_total, + entropy_modified, + energy_total_modified, + lake_at_rest_error)) + @test_nowarn print(analysis_callback) + @test_nowarn display(analysis_callback) + end + + @testset "RelaxationCallback" begin + relaxation_callback = RelaxationCallback(invariant = entropy) + @test_nowarn print(relaxation_callback) + @test_nowarn display(relaxation_callback) + end + + @testset "util" begin + @test_nowarn get_examples() + @test_nowarn DispersiveShallowWater.path_create_figures() + @test_nowarn trixi_include(default_example(), tspan = (0.0, 0.1)) + + accuracy_orders = [2, 4, 6] + for accuracy_order in accuracy_orders + eoc_mean_values, _ = convergence_test(default_example(), 2, initial_N = 128, + tspan = (0.0, 1.0), + accuracy_order = accuracy_order) + show(eoc_mean_values) + @test isapprox(eoc_mean_values[:l2][1], accuracy_order, atol = 0.5) + @test isapprox(eoc_mean_values[:linf][2], accuracy_order, atol = 0.5) + @test isapprox(eoc_mean_values[:l2][1], accuracy_order, atol = 0.5) + @test isapprox(eoc_mean_values[:linf][2], accuracy_order, atol = 0.5) + end end end diff --git a/test/test_visualization.jl b/test/test_visualization.jl new file mode 100644 index 00000000..7ffaac8e --- /dev/null +++ b/test/test_visualization.jl @@ -0,0 +1,21 @@ +module TestUnit + +using Test +using DispersiveShallowWater +using Plots + +@testset "Unit tests" begin + @testset "Visualization" begin + trixi_include(@__MODULE__, default_example(), tspan = (0.0, 1.0)) + @test_nowarn plot(semi => sol) + @test_nowarn plot!(semi => sol, plot_initial = true) + @test_nowarn plot(semi, sol, conversion = prim2cons, plot_bathymetry = false) + @test_nowarn plot(semi => sol, 0.0) + @test_nowarn plot(semi, sol, 0.0, conversion = prim2cons) + @test_nowarn plot(analysis_callback) + @test_nowarn plot(analysis_callback, what = (:errors,)) + @test_nowarn plot(analysis_callback, what = (:integrals, :errors)) + end +end + +end # module diff --git a/visualization/create_figures.jl b/visualization/create_figures.jl index f2165da6..873cc933 100644 --- a/visualization/create_figures.jl +++ b/visualization/create_figures.jl @@ -286,7 +286,7 @@ function fig_7() c = 5 / 2 * sqrt(g * D) xmin = -35.0 xmax = 35.0 - tspan = (0.0, 20 * (xmax - xmin) / c) + tspan = (0.0, 15 * (xmax - xmin) / c) N = 512 accuracy_order = 8