Manufactured solution for Svärd-Kalisch equations with variable bathy…

…metry (#84) * test with variable bathymetry for manufactured solution * WIP: different tries for manufactured solution of Svärd-Kalisch equations with variable bathymetry * try to fix manufactured solution with variable bathymetry * finally get manufactured solution for Svärd-Kalisch equations with variable bathymetry to work * add test * format * Update src/equations/svaerd_kalisch_1d.jl Co-authored-by: Hendrik Ranocha <[email protected]> * lower tolerance to let CI pass * fix location of tolerance * increase interval of AnalysisCallback --------- Co-authored-by: Hendrik Ranocha <[email protected]>
JoshuaLampert · Feb 2, 2024 · ccb42cb · ccb42cb
1 parent 8946fea
commit ccb42cb
Show file tree

Hide file tree

Showing 5 changed files with 62 additions and 32 deletions.
diff --git a/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_manufactured.jl b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_manufactured.jl
@@ -33,12 +33,12 @@ semi = Semidiscretization(mesh, equations, initial_condition, solver,
 tspan = (0.0, 1.0)
 ode = semidiscretize(semi, tspan)
 summary_callback = SummaryCallback()
-analysis_callback = AnalysisCallback(semi; interval = 10,
+analysis_callback = AnalysisCallback(semi; interval = 1000,
                                      extra_analysis_errors = (:conservation_error,),
                                      extra_analysis_integrals = (waterheight_total,
                                                                  velocity, entropy))
 callbacks = CallbackSet(analysis_callback, summary_callback)
 
 saveat = range(tspan..., length = 100)
-sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7,
+sol = solve(ode, Tsit5(), abstol = 1e-12, reltol = 1e-12,
             save_everystep = false, callback = callbacks, saveat = saveat)
diff --git a/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_well_balanced.jl b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_well_balanced.jl
@@ -51,7 +51,7 @@ semi = Semidiscretization(mesh, equations, initial_condition, solver,
 tspan = (0.0, 10.0)
 ode = semidiscretize(semi, tspan)
 summary_callback = SummaryCallback()
-analysis_callback = AnalysisCallback(semi; interval = 10,
+analysis_callback = AnalysisCallback(semi; interval = 100,
                                      extra_analysis_errors = (:conservation_error,),
                                      extra_analysis_integrals = (waterheight_total,
                                                                  momentum, entropy,

diff --git a/src/equations/bbm_bbm_variable_bathymetry_1d.jl b/src/equations/bbm_bbm_variable_bathymetry_1d.jl
@@ -67,7 +67,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 = 5.0 + 2.0 * cospi(2 * x)
+    D = 5 + 2 * cospi(2 * x)
     return SVector(eta, v, D)
 end
 
@@ -97,7 +97,7 @@ function source_terms_manufactured(q, x, t, equations::BBMBBMVariableEquations1D
            4 * (a4 + 2 * pi * a3) * (16 * sinpi(x)^4 - 26 * sinpi(x)^2 + 7)) * exp(t / 2) /
           3 - exp(t / 2) * a4 / 2 - pi * exp(t / 2) * a3 - pi * exp(t) * a5
 
-    return SVector(dq1, dq2, 0.0)
+    return SVector(dq1, dq2, zero(dq1))
 end
 
 """

diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl
@@ -94,7 +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 = 3.0
+    D = 5 + 2 * cospi(2 * x)
     return SVector(eta, v, D)
 end
 
@@ -105,31 +105,61 @@ 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 is constant, thus simply take the first entry
     eta0 = equations.eta0
     alpha = equations.alpha
     beta = equations.beta
     gamma = equations.gamma
-    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)
+    a1 = sinpi(2 * x)
+    a2 = cospi(2 * x)
+    a3 = sinpi(-t + 2 * x)
+    a4 = cospi(-t + 2 * x)
+    a5 = sinpi(t - 2 * x)
+    a6 = cospi(t - 2 * x)
+    a7 = sinpi(-4 * t + 2 * x)
+    a8 = exp(t / 2)
+    a9 = exp(t) * cospi(-4 * t + 2 * x)
+    a10 = eta0 + 2.0 * a2 + 5.0
+    a11 = sqrt(g * a10)
+    a12 = 0.2 * eta0 + 0.4 * a2 + 1
+    a13 = alpha * a11 * a12^2
+    a14 = sqrt(a13)
+    a15 = -1.0 * pi * a13 * a1 / a10 - 0.8 * pi * alpha * a11 * a12 * a1
+    a16 = -20.0 * pi^2 * a14 * a9 - 10.0 * pi * a14 * a15 * exp(t) * a7 / (a13)
+    a17 = -2 * pi * exp(t) * a7 - 4.0 * pi * a1
+    a18 = a9 + 2.0 * a2 + 5.0
+    a19 = a17 * a8 * a3 + 2 * pi * a18 * a8 * a4
+    a20 = a14 * (40.0 * pi^3 * a14 * exp(t) * a7 - 40.0 * pi^2 * a14 * a15 * a9 / (a13) -
+           20.0 * pi^2 * a14 * a15 * exp(t) * a1 * a7 / (a13 * a10) -
+           16.0 * pi^2 * a14 * a15 * exp(t) * a1 * a7 / (alpha * a11 * a12^3) -
+           10.0 * pi * a14 *
+           (-2.0 * pi^2 * a13 * a2 / a10 - 1.6 * pi^2 * alpha * a11 * a12 * a2 +
+            3.2 * pi^2 * alpha * a11 * a12 * a1^2 / a10 + 0.56 * pi^2 * alpha * a11 * a1^2) *
+           exp(t) * a7 / (a13) -
+           10.0 * pi * a14 * a15^2 * exp(t) * a7 / (alpha^2 * g * a12^4 * a10))
+
+    dq1 = -5.0 * a20 + a19 + 4 * pi * exp(t) * a7 + a9 - 5.0 * a14 * a16 * a15 / (a13)
+
+    dq2 = -25.0 * beta * (-2 * pi^2 * a8 * a3 + 4 * pi^3 * a8 * a4) * a12^2 * a10 +
+          100.0 * pi * beta * (2 * pi^2 * a8 * a3 + pi * a8 * a4) * a12^2 * a1 +
+          40.0 * pi * beta * (2 * pi^2 * a8 * a3 + pi * a8 * a4) * a12 * a10 * a1 -
+          2 * pi * g * a18 * exp(t) * a7 +
+          100.0 * pi^3 * gamma * a11 * a12^2 * a10 * a8 * a4 -
+          300.0 * pi^3 * gamma * a11 * a12^2 * a8 * a1 * a3 -
+          80.0 * pi^3 * gamma * a11 * a12 * a10 * a8 * a1 * a3 -
+          pi^3 * gamma * a11 *
+          (-50.0 * (3.2 * a12 * a2 - 1.28 * a1^2) * a10 * a6 -
+           50.0 * (4.0 * a2 / a10 + 0.16 * a1^2 / a12^2) * a12^2 * a10 * a6 -
+           200.0 * a12^2 * a10 * a6 - 1200.0 * a12^2 * a1 * a5 - 400.0 * a12^2 * a2 * a6 +
+           800.0 * a12^2 * a1^2 * a6 / a10 - 320.0 * a12 * a10 * a1 * a5 +
+           960.0 * a12 * a1^2 * a6) * a8 / 2 - 10.0 * pi * a14 * a16 * a8 * a4 -
+          2.5 * a20 * a8 * a3 + (5.0 * a20 + 5.0 * a14 * a16 * a15 / (a13)) * a8 * a3 / 2 +
+          (a8 * a3 / 2 - pi * a8 * a4) * a18 + a17 * exp(t) * a3^2 / 2 -
+          (a19) * a8 * a3 / 2 + 3 * pi * a18 * exp(t) * a3 * a4 -
+          2.5 * a14 * a16 * a15 * a8 * a3 / (a13)
+
+    return SVector(dq1, dq2, zero(dq1))
 end
-#
+
 function create_cache(mesh,
                       equations::SvaerdKalischEquations1D,
                       solver,

diff --git a/test/test_svaerd_kalisch_1d.jl b/test/test_svaerd_kalisch_1d.jl
@@ -12,12 +12,12 @@ EXAMPLES_DIR = joinpath(examples_dir(), "svaerd_kalisch_1d")
         @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,
-                            atol=1e-4) # in order to make CI pass
+                            l2=[3.3755102050606554e-6 2.320896961343125e-7 0.0],
+                            linf=[4.908886917176503e-6 3.888671399332466e-7 0.0],
+                            cons_error=[2.42861286636753e-16 1.9224170696150768e-7 0.0],
+                            change_waterheight=-2.42861286636753e-16,
+                            change_entropy=0.1868146724821993,
+                            atol=1e-9) # in order to make CI pass
     end
 
     @trixi_testset "svaerd_kalisch_1d_dingemans" begin