From 672cd3ec698e2b771ed8546efb71790ef6ea3437 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Mon, 2 Oct 2023 03:51:31 +0200 Subject: [PATCH] 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