add more unit tests and CG elixir

JoshuaLampert · Oct 2, 2023 · 672cd3e · 672cd3e
1 parent 45da30b
commit 672cd3e
Show file tree

Hide file tree

Showing 3 changed files with 108 additions and 6 deletions.
diff --git a/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans_cg.jl 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
@@ -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
@@ -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,14 +124,15 @@ 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
         equations = @test_nowarn SvaerdKalischEquations1D(gravity_constant = 9.81,
                                                           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