add support of source terms for Svärd-Kalisch equations with constant…

… bathymetry
JoshuaLampert · Nov 19, 2023 · ffd5fb5 · ffd5fb5
1 parent 238f704
commit ffd5fb5
Show file tree

Hide file tree

Showing 3 changed files with 34 additions and 17 deletions.
diff --git a/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_manufactured.jl b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_manufactured.jl
@@ -2,9 +2,9 @@ using OrdinaryDiffEq
 using DispersiveShallowWater
 
 ###############################################################################
-# Semidiscretization of the BBM-BBM equations
+# Semidiscretization of the Svärd-Kalisch equations
 
-equations = SvaerdKalischEquations1D(gravity_constant = 1.0, eta0 = 2.0,
+equations = SvaerdKalischEquations1D(gravity_constant = 1.0, eta0 = 0.0,
                                      alpha = 0.0004040404040404049,
                                      beta = 0.49292929292929294,
                                      gamma = 0.15707070707070708)
@@ -16,7 +16,7 @@ boundary_conditions = boundary_condition_periodic
 # create homogeneous mesh
 coordinates_min = 0.0
 coordinates_max = 1.0
-N = 512
+N = 128
 mesh = Mesh1D(coordinates_min, coordinates_max, N)
 
 # create solver with periodic SBP operators of accuracy order 4
@@ -40,4 +40,4 @@ callbacks = CallbackSet(analysis_callback)
 
 saveat = range(tspan..., length = 100)
 sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7,
-            save_everystep = true, callback = callbacks, saveat = saveat)
+            save_everystep = false, callback = callbacks, saveat = saveat)
diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl
@@ -94,8 +94,7 @@ function initial_condition_manufactured(x, t,
                                         mesh)
     eta = exp(t) * cospi(2 * (x - 2 * t))
     v = exp(t / 2) * sinpi(2 * (x - t / 2))
-#     D = cospi(2 * x)
-    D = -2.0
+    D = 3.0
     return SVector(eta, v, D)
 end
 
@@ -106,21 +105,28 @@ A smooth manufactured solution in combination with [`initial_condition_manufactu
 """
 function source_terms_manufactured(q, x, t, equations::SvaerdKalischEquations1D)
     g = equations.gravity
-    D = q[3, 1]
+    D = q[3, 1] # D is constant, thus simply take the first entry
     eta0 = equations.eta0
     alpha = equations.alpha
     beta = equations.beta
     gamma = equations.gamma
-    a1 = cospi(2 * x)
-    a2 = sinpi(2 * x)
-    a3 = cospi(t - 2 * x)
-    a4 = sinpi(t - 2 * x)
-    a6 = sinpi(4 * t - 2 * x)
-    a7 = cospi(4 * t - 2 * x)
-#     dq1 = -10*pi^3*alpha^2*g*(eta0 + a1)^3*(2*(eta0 + a1)*a7 + 5*a2*a6)*exp(t)*a2 + 2*pi^3*alpha^2*g*(eta0 + a1)^3*(4*(eta0 + a1)^2*a6 - 20*(eta0 + a1)*a2*a7 + 10*(eta0 + a1)*a6*a1 - 15*a2^2*a6)*exp(t) - 2*pi*(exp(t)*a6 - a2)*exp(t/2)*a4 + 2*pi*(exp(t)*a7 + a1)*exp(t/2)*a3 - 4*pi*exp(t)*a6 + exp(t)*a7
-#     dq2 = 4*pi^3*alpha^2*g*(eta0 + a1)^4*(2*(eta0 + a1)*a7 + 5*a2*a6)*exp(3*t/2)*a3 + 10*pi^3*alpha^2*g*(eta0 + a1)^3*(2*(eta0 + a1)*a7 + 5*a2*a6)*exp(3*t/2)*a2*a4 - 2*pi^3*alpha^2*g*(eta0 + a1)^3*(4*(eta0 + a1)^2*a6 - 20*(eta0 + a1)*a2*a7 + 10*(eta0 + a1)*a6*a1 - 15*a2^2*a6)*exp(3*t/2)*a4 - 2*pi^2*beta*(eta0 + a1)^2*((eta0 + a1)*a4 + 2*pi*(eta0 + a1)*a3 + 6*pi*a2*a4 - 3*a2*a3)*exp(t/2) + 2*pi*g*(exp(t)*a7 + a1)*exp(t)*a6 + 4.0*pi^3*gamma*sqrt(g*(eta0 + a1))*(eta0 + a1)^3*exp(t/2)*a3 + 14.0*pi^3*gamma*sqrt(g*(eta0 + a1))*(eta0 + a1)^2*exp(t/2)*a2*a4 + pi^3*gamma*sqrt(g*(eta0 + a1))*(eta0 + a1)*(4*(eta0 + a1)^2*a3 + 28*(eta0 + a1)*a2*a4 + 12*((eta0 + a1)*a1 - 2*a2^2)*a3 + (2*(eta0 + a1)*a1 + a2^2)*a3 - 12*a2^2*a3)*exp(t/2) + (4*pi*a6 - a7)*exp(3*t/2)*a4 + 2*pi*(exp(t)*a6 - a2)*exp(t)*a4^2 - (exp(t)*a7 + a1)*exp(t/2)*a4/2 - pi*(exp(t)*a7 + a1)*exp(t/2)*a3 - 4*pi*(exp(t)*a7 + a1)*exp(t)*a4*a3 - (10*pi^3*alpha^2*g*(eta0 + a1)^3*(2*(eta0 + a1)*a7 + 5*a2*a6)*exp(t)*a2 - 2*pi^3*alpha^2*g*(eta0 + a1)^3*(4*(eta0 + a1)^2*a6 - 20*(eta0 + a1)*a2*a7 + 10*(eta0 + a1)*a6*a1 - 15*a2^2*a6)*exp(t) + 2*pi*(exp(t)*a6 - a2)*exp(t/2)*a4 - 2*pi*(exp(t)*a7 + a1)*exp(t/2)*a3 + 4*pi*exp(t)*a6 - exp(t)*a7)*exp(t/2)*a4
-    dq1 = 8*pi^3*alpha^2*g*(D + eta0)^5*exp(t)*sin(pi*(4*t - 2*x)) + 2*pi*(D + exp(t)*cos(pi*(4*t - 2*x)))*exp(t/2)*cos(pi*(t - 2*x)) - 2*pi*exp(3*t/2)*sin(pi*(t - 2*x))*sin(pi*(4*t - 2*x)) - 4*pi*exp(t)*sin(pi*(4*t - 2*x)) + exp(t)*cos(pi*(4*t - 2*x))
-    dq2 = 2*pi*D*g*exp(t)*sin(4*pi*t - 2*pi*x) - D*exp(t/2)*sin(pi*t - 2*pi*x)/2 - pi*D*exp(t/2)*cos(pi*t - 2*pi*x) - 2*pi*D*exp(t)*sin(pi*t - 2*pi*x)*cos(pi*t - 2*pi*x) + 8*pi^3*alpha^2*g*(D + eta0)^5*exp(3*t/2)*cos(pi*t - 2*pi*x)*cos(4*pi*t - 2*pi*x) - 2*pi^2*beta*(D + eta0)^3*exp(t/2)*sin(pi*t - 2*pi*x) - 4*pi^3*beta*(D + eta0)^3*exp(t/2)*cos(pi*t - 2*pi*x) + 2*pi*g*exp(2*t)*sin(4*pi*t - 2*pi*x)*cos(4*pi*t - 2*pi*x) + 8.0*pi^3*gamma*(D + eta0)^3*sqrt(D*g + eta0*g)*exp(t/2)*cos(pi*t - 2*pi*x) - exp(3*t/2)*sin(pi*t - 2*pi*x)*cos(4*pi*t - 2*pi*x)/2 - pi*exp(3*t/2)*cos(pi*t - 2*pi*x)*cos(4*pi*t - 2*pi*x) - 2*pi*exp(2*t)*sin(pi*t - 2*pi*x)*cos(pi*t - 2*pi*x)*cos(4*pi*t - 2*pi*x)
+    a1 = cospi(t - 2 * x)
+    a2 = sinpi(t - 2 * x)
+    a3 = cospi(4 * t - 2 * x)
+    a4 = sinpi(4 * t - 2 * x)
+
+    dq1 = 8 * pi^3 * alpha * sqrt(g * (D + eta0)) * (D + eta0)^2 * exp(t) * a4 +
+          2 * pi * (D + exp(t) * a3) * exp(t / 2) * a1 - 2 * pi * exp(3 * t / 2) * a2 * a4 -
+          4 * pi * exp(t) * a4 + exp(t) * a3
+    dq2 = 2 * pi * D * g * exp(t) * a4 - D * exp(t / 2) * a2 / 2 -
+          pi * D * exp(t / 2) * a1 - 2 * pi * D * exp(t) * a2 * a1 +
+          8 * pi^3 * alpha * (D + eta0)^2 * sqrt(D * g + eta0 * g) * exp(3 * t / 2) * a1 *
+          a3 - 2 * pi^2 * beta * (D + eta0)^3 * exp(t / 2) * a2 -
+          4 * pi^3 * beta * (D + eta0)^3 * exp(t / 2) * a1 +
+          2 * pi * g * exp(2 * t) * a4 * a3 +
+          8.0 * pi^3 * gamma * (D + eta0)^3 * sqrt(D * g + eta0 * g) * exp(t / 2) * a1 -
+          exp(3 * t / 2) * a2 * a3 / 2 - pi * exp(3 * t / 2) * a1 * a3 -
+          2 * pi * exp(2 * t) * a2 * a1 * a3
     return SVector(dq1, dq2, 0.0)
 end
 #

diff --git a/test/test_svaerd_kalisch_1d.jl b/test/test_svaerd_kalisch_1d.jl
@@ -8,6 +8,17 @@ include("test_util.jl")
 EXAMPLES_DIR = joinpath(examples_dir(), "svaerd_kalisch_1d")
 
 @testset "SvaerdKalisch1D" begin
+    @trixi_testset "svaerd_kalisch_1d_manufactured" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                     "svaerd_kalisch_1d_manufactured.jl"),
+                            tspan=(0.0, 0.1),
+                            l2=[7.467887060263923e-5 2.7796353838948894e-8 0.0],
+                            linf=[0.0001613144395267163 4.344495230235168e-8 0.0],
+                            cons_error=[2.3635607360183997e-16 8.084235123776567e-10 0.0],
+                            change_waterheight=-2.3635607360183997e-16,
+                            change_entropy=0.1342289500320556)
+    end
+
     @trixi_testset "svaerd_kalisch_1d_dingemans" begin
         @test_trixi_include(joinpath(EXAMPLES_DIR,
                                      "svaerd_kalisch_1d_dingemans.jl"),