From 672cd3ec698e2b771ed8546efb71790ef6ea3437 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Mon, 2 Oct 2023 03:51:31 +0200 Subject: [PATCH 1/5] add more unit tests and CG elixir --- .../svaerd_kalisch_1d_dingemans_cg.jl | 47 ++++++++++++++++ test/test_svaerd_kalisch_1d.jl | 12 ++++ test/test_unit.jl | 55 +++++++++++++++++-- 3 files changed, 108 insertions(+), 6 deletions(-) create mode 100644 examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans_cg.jl 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/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..cc25cc31 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -7,31 +7,36 @@ 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 +59,47 @@ 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) q = [42.0, 2.0] @test prim2prim(q, equations) == q @test isapprox(cons2prim(prim2cons(q, equations), equations), q) @@ -67,11 +109,12 @@ 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) q = [42.0, 2.0, 2.0] @test prim2prim(q, equations) == q @test isapprox(cons2prim(prim2cons(q, equations), equations), q) @@ -81,7 +124,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 +131,8 @@ using SparseArrays: sparse, SparseMatrixCSC alpha = 0.0004040404040404049, beta = 0.49292929292929294, gamma = 0.15707070707070708) + @test_nowarn print(equations) + @test_nowarn display(equations) q = [42.0, 2.0, 2.0] @test prim2prim(q, equations) == q @test isapprox(cons2prim(prim2cons(q, equations), equations), q) @@ -98,7 +142,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 end From 63f3364cba53550b185ff86319500573f8fee181 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Mon, 2 Oct 2023 11:11:31 +0200 Subject: [PATCH 2/5] format --- test/test_unit.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/test_unit.jl b/test/test_unit.jl index cc25cc31..abe56c21 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -37,7 +37,6 @@ using SparseArrays: sparse, SparseMatrixCSC @test grid(solver) == grid(solver.D1) == grid(solver.D2) @test real(solver) == Float64 - 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) From dce8a3e857fc0a1057be9bdd3415a026e9c26674 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Mon, 2 Oct 2023 11:47:21 +0200 Subject: [PATCH 3/5] more unit tests --- src/callbacks_step/analysis.jl | 2 +- test/test_unit.jl | 88 ++++++++++++++++++++++++++++++++++ 2 files changed, 89 insertions(+), 1 deletion(-) 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/test_unit.jl b/test/test_unit.jl index abe56c21..1ca46be8 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -99,6 +99,20 @@ using SparseArrays: sparse, SparseMatrixCSC 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) @@ -114,6 +128,20 @@ using SparseArrays: sparse, SparseMatrixCSC 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) @@ -132,6 +160,22 @@ using SparseArrays: sparse, SparseMatrixCSC 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) @@ -142,6 +186,50 @@ using SparseArrays: sparse, SparseMatrixCSC @test discharge(q, equations) == 88.0 @test isapprox(energy_total(q, equations), 8740.42) 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 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 end # module From 22e1b068164f0f42d35fb0c618c679c5fa3645ab Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Mon, 2 Oct 2023 12:23:04 +0200 Subject: [PATCH 4/5] test visualization --- test/Project.toml | 2 ++ test/runtests.jl | 1 + test/test_unit.jl | 1 + test/test_visualization.jl | 21 +++++++++++++++++++++ visualization/create_figures.jl | 2 +- 5 files changed, 26 insertions(+), 1 deletion(-) create mode 100644 test/test_visualization.jl 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_unit.jl b/test/test_unit.jl index 1ca46be8..ebe64afc 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -216,6 +216,7 @@ using SparseArrays: sparse, SparseMatrixCSC @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] diff --git a/test/test_visualization.jl b/test/test_visualization.jl new file mode 100644 index 00000000..7a1a3805 --- /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 From 79014e6fbf9f63b11cadc19c537af0c5751bfd23 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Mon, 2 Oct 2023 12:23:52 +0200 Subject: [PATCH 5/5] format --- test/test_visualization.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_visualization.jl b/test/test_visualization.jl index 7a1a3805..7ffaac8e 100644 --- a/test/test_visualization.jl +++ b/test/test_visualization.jl @@ -14,7 +14,7 @@ using Plots @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,)) + @test_nowarn plot(analysis_callback, what = (:integrals, :errors)) end end