Fix ShallowWaterExnerEquations1D for negative / zero velocities (#101)

* compare against absolute values for the velocity

* add symmetric dam break test case

* Setup type conversion for sediment models and add type tests

* apply formatter

* update test values due to updated roe average for h_s

* pass correct sediment height to eigenvalue computation

Author: patrickersing

examples/tree_1d_dgsem/elixir_shallowwater_exner_channel.jl

@@ -12,7 +12,7 @@ equations = ShallowWaterExnerEquations1D(gravity = 9.81, rho_f = 1.0, rho_s = 1.
                                          porosity = 0.4,
                                          sediment_model = GrassModel(A_g = 0.01))
 
-# Initial condition for for a channel flow problem over a sediment hump.
+# Initial condition for a channel flow problem over a sediment hump.
 # An asymptotic solution based on the method of characteristics was derived under a rigid-lid
 # approximation in chapter 3.5.1 of the thesis:
 # - Justin Hudson (2001)

examples/tree_1d_dgsem/elixir_shallowwater_exner_dam_break_symmetric.jl

@@ -0,0 +1,99 @@
+
+using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the shallow water exner equations for a subcritical symmetric dam break problem
+
+equations = ShallowWaterExnerEquations1D(gravity = 9.81, rho_f = 0.3, rho_s = 1.0,
+                                         porosity = 0.4,
+                                         sediment_model = GrassModel(A_g = 0.01))
+
+# Initial condition for a synthetic test case of a symmetric dam break problem.
+# The setup is taken from the paper:
+# - S. Martínez-Aranda, J. Murillo, P. García-Navarro (2021)
+#   "Comparison of new efficient 2D models for the simulation of bedload transport using the augmented
+#    Roe approach"
+#   [DOI: 10.1016/j.advwatres.2021.103931](https://doi.org/10.1016/j.advwatres.2021.103931)
+function initial_condition_dam_break_symmetric(x, t,
+                                               equations::ShallowWaterExnerEquations1D)
+    # Setup initial perturbation of the water height
+    if -0.5 <= x[1] <= 0.5
+        h = 1.0
+    else
+        h = 0.2
+    end
+
+    # Set constant values for the sediment height and zero momentum
+    hv = 0.0
+    h_b = 1.0
+
+    return SVector(h, hv, h_b)
+end
+
+initial_condition = initial_condition_dam_break_symmetric
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+surface_flux = (FluxPlusDissipation(flux_ersing_etal, dissipation_roe),
+                flux_nonconservative_ersing_etal)
+
+basis = LobattoLegendreBasis(3)
+
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max = 0.5,
+                                         alpha_min = 0.001,
+                                         alpha_smooth = true,
+                                         variable = water_sediment_height)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+###############################################################################
+# Get the TreeMesh and setup a periodic mesh
+
+coordinates_min = -4.0
+coordinates_max = 4.0
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 7,
+                n_cells_max = 10_000,
+                periodicity = true)
+
+# Create the semi discretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solver
+
+tspan = (0.0, 0.6)
+ode = semidiscretize(semi, tspan)
+
+###############################################################################
+# Callbacks
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 1.0)
+
+save_solution = SaveSolutionCallback(dt = 0.1,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback,
+                        stepsize_callback, save_solution)
+
+###############################################################################
+# run the simulation
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false);
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            ode_default_options()..., callback = callbacks);

src/equations/shallow_water_exner_1d.jl

@@ -80,8 +80,9 @@ struct GrassModel{RealT} <: SedimentModel{RealT}
     m_g::RealT
 end
 
-function GrassModel(; A_g, m_g = 3.0)
-    return GrassModel(A_g, m_g)
+function GrassModel(; A_g, m_g = 3)
+    RealT = promote_type(typeof(A_g), typeof(m_g))
+    return GrassModel(RealT(A_g), RealT(m_g))
 end
 
 @doc raw"""
@@ -97,7 +98,9 @@ An overview of different formulations to compute the sediment discharge can be f
   [DOI:10.1016/j.compfluid.2007.07.017](https://doi.org/10.1016/j.compfluid.2007.07.017)
 """
 function MeyerPeterMueller(; theta_c, d_s)
-    return ShieldsStressModel(0.0, 1.5, 0.0, 8.0, 1.0, 0.0, theta_c, d_s)
+    RealT = promote_type(typeof(theta_c), typeof(d_s))
+    return ShieldsStressModel(RealT(0.0), RealT(1.5), RealT(0.0), RealT(8.0),
+                              RealT(1.0), RealT(0.0), RealT(theta_c), RealT(d_s))
 end
 
 @doc raw"""
@@ -228,9 +231,9 @@ A smooth initial condition used for convergence tests in combination with
                                                           equations::ShallowWaterExnerEquations1D)
     ω = sqrt(2) * pi
 
-    h = 2.0 + cos(ω * x[1]) * cos(ω * t)
-    v = 0.5
-    h_b = 2.0 + sin(ω * x[1]) * cos(ω * t)
+    h = 2 + cos(ω * x[1]) * cos(ω * t)
+    v = 0.5f0
+    h_b = 2 + sin(ω * x[1]) * cos(ω * t)
 
     return SVector(h, h * v, h_b)
 end
@@ -256,20 +259,21 @@ equations = ShallowWaterExnerEquations1D(gravity = 10.0, rho_f = 0.5,
                                                                                                                     T,
                                                                                                                     S
                                                                                                                     }
-    ω = sqrt(2.0) * pi
+    ω = sqrt(2) * pi
     A_g = equations.sediment_model.A_g
 
-    h = -cos(x[1] * ω) * sin(t * ω) * ω - 0.5 * sin(x[1] * ω) * cos(t * ω) * ω
-    hv = -0.5 * cos(x[1] * ω) * sin(t * ω) * ω - 0.25 * sin(x[1] * ω) * cos(t * ω) * ω +
-         10.0 * A_g *
-         (cos(x[1] * ω) * cos(t * ω) * ω - 0.5 * sin(x[1] * ω) * cos(t * ω) * ω) +
-         10.0 * (2.0 + cos(x[1] * ω) * cos(t * ω)) *
+    h = -cos(x[1] * ω) * sin(t * ω) * ω - 0.5f0 * sin(x[1] * ω) * cos(t * ω) * ω
+    hv = -0.5f0 * cos(x[1] * ω) * sin(t * ω) * ω -
+         0.25f0 * sin(x[1] * ω) * cos(t * ω) * ω +
+         10 * A_g *
+         (cos(x[1] * ω) * cos(t * ω) * ω - 0.5f0 * sin(x[1] * ω) * cos(t * ω) * ω) +
+         10 * (2 + cos(x[1] * ω) * cos(t * ω)) *
          (cos(x[1] * ω) * cos(t * ω) * ω - sin(x[1] * ω) * cos(t * ω) * ω)
     h_b = -sin(x[1] * ω) * sin(t * ω) * ω
     return SVector(h, hv, h_b)
 end
 
-"""
+""" 
     source_terms_convergence_test(u, x, t, equations::ShallowWaterExnerEquations1D{T, S, ShieldsStressModel{T}}) where {T, S}
 
 Source terms used for convergence tests in combination with [`Trixi.initial_condition_convergence_test`](@ref)
@@ -291,27 +295,27 @@ equations = ShallowWaterExnerEquations1D(gravity = 10.0, rho_f = 0.5,
                                                                                                                             T,
                                                                                                                             S
                                                                                                                             }
-    ω = sqrt(2.0) * pi
+    ω = sqrt(2) * pi
     (; gravity, porosity_inv, rho_f, rho_s, r) = equations
 
     n = equations.friction.n
 
     # Constant expression from the MPM model
-    c = sqrt(gravity * (1 / r - 1)) * 8.0 * porosity_inv *
+    c = sqrt(gravity * (1 / r - 1)) * 8 * porosity_inv *
         (rho_f / (rho_s - rho_f))^(3 / 2) * n^3
 
-    h = -cos(x[1] * ω) * sin(t * ω) * ω - 0.5 * sin(x[1] * ω) * cos(t * ω) * ω
+    h = -cos(x[1] * ω) * sin(t * ω) * ω - 0.5f0 * sin(x[1] * ω) * cos(t * ω) * ω
 
-    hv = ((5.0 * c *
-           (cos(x[1] * ω) * cos(t * ω) * ω - 0.5 * sin(x[1] * ω) * cos(t * ω) * ω)) /
-          ((2.0 + cos(x[1] * ω) * cos(t * ω))^0.5) -
-          0.5 * cos(x[1] * ω) * sin(t * ω) * ω -
-          0.25 * sin(x[1] * ω) * cos(t * ω) * ω +
-          10.0 * (2.0 + cos(x[1] * ω) * cos(t * ω)) *
+    hv = ((5 * c *
+           (cos(x[1] * ω) * cos(t * ω) * ω - 0.5f0 * sin(x[1] * ω) * cos(t * ω) * ω)) /
+          ((2 + cos(x[1] * ω) * cos(t * ω))^0.5) -
+          0.5f0 * cos(x[1] * ω) * sin(t * ω) * ω -
+          0.25f0 * sin(x[1] * ω) * cos(t * ω) * ω +
+          10 * (2 + cos(x[1] * ω) * cos(t * ω)) *
           (cos(x[1] * ω) * cos(t * ω) * ω - sin(x[1] * ω) * cos(t * ω) * ω))
 
-    h_b = ((0.5 * ((0.125 * c) / (2.0 + cos(x[1] * ω) * cos(t * ω))) * sin(x[1] * ω) *
-            cos(t * ω) * ω) / ((2.0 + cos(x[1] * ω) * cos(t * ω))^0.5) -
+    h_b = ((0.5f0 * ((0.125f0 * c) / (2 + cos(x[1] * ω) * cos(t * ω))) * sin(x[1] * ω) *
+            cos(t * ω) * ω) / ((2 + cos(x[1] * ω) * cos(t * ω))^0.5) -
            sin(x[1] * ω) * sin(t * ω) * ω)
 
     return SVector(h, hv, h_b)
@@ -367,8 +371,8 @@ scheme that is entropy conservative and well-balanced.
     h_b_jump = h_b_rr - h_b_ll
 
     # Workaround to avoid division by zero, when computing the effective sediment height
-    if velocity(u_ll, equations) < eps()
-        h_s_ll = 0.0
+    if abs(velocity(u_ll, equations)) < eps(typeof(h_ll))
+        h_s_ll = 0
     else
         h_s_ll = q_s(u_ll, equations) / velocity(u_ll, equations)
     end
@@ -405,12 +409,12 @@ To obtain an entropy stable formulation the `surface_flux` can be set as
     v_rr = velocity(u_rr, equations)
 
     # Average each factor of products in flux
-    v_avg = 0.5 * (v_ll + v_rr)
+    v_avg = 0.5f0 * (v_ll + v_rr)
 
     # Calculate fluxes depending on orientation
-    f1 = 0.5 * (h_v_ll + h_v_rr)
+    f1 = 0.5f0 * (h_v_ll + h_v_rr)
     f2 = f1 * v_avg
-    f3 = 0.5 * (q_s(u_ll, equations) + q_s(u_rr, equations))
+    f3 = 0.5f0 * (q_s(u_ll, equations) + q_s(u_rr, equations))
 
     return SVector(f1, f2, f3)
 end
@@ -432,21 +436,24 @@ for the sediment discharge `q_s`.
     v_rr = velocity(u_rr, equations)
 
     # Compute approximate Roe averages.
-    # The actual Roe average for the sediment discharge `q_s` would depend on the sediment and
-    # friction model and can be difficult to compute analytically.
+    # The actual Roe average for the sediment height `h_b` depends on the sediment and
+    # friction model and an explicit formula is not always available.
     # Therefore we only use an approximation here.
-    h_avg = 0.5 * (u_ll[1] + u_rr[1])
+    h_avg = 0.5f0 * (u_ll[1] + u_rr[1])
     v_avg = (sqrt(u_ll[1]) * v_ll + sqrt(u_rr[1]) * v_rr) /
             (sqrt(u_ll[1]) + sqrt(u_rr[1]))
+    h_b_avg = 0.5f0 * (u_ll[3] + u_rr[3])
+
     # Workaround to avoid division by zero, when computing the effective sediment height
-    if v_avg < eps()
-        h_s_avg = 0.0
+    if abs(v_avg) < eps(typeof(h_avg))
+        h_s_avg = 0
     else
-        h_s_avg = 0.5 * (q_s(u_ll, equations) / v_ll + q_s(u_rr, equations) / v_rr)
+        h_s_avg = (q_s(SVector(h_avg, h_avg * v_avg, h_b_avg), equations) /
+                   v_avg)
     end
 
     # Compute the eigenvalues using Cardano's formula
-    λ1, λ2, λ3 = eigvals_cardano(SVector(h_avg, h_avg * v_avg, h_s_avg), equations)
+    λ1, λ2, λ3 = eigvals_cardano(SVector(h_avg, h_avg * v_avg, h_b_avg), equations)
 
     # Precompute some common expressions
     c1 = g * (h_avg + h_s_avg)
@@ -470,7 +477,7 @@ for the sediment discharge `q_s`.
     Λ_abs = @SMatrix [abs(λ1) z z; z abs(λ2) z; z z abs(λ3)]
 
     # Compute dissipation
-    diss = SVector(-0.5 * R * Λ_abs * R_inv * (u_rr - u_ll))
+    diss = SVector(-0.5f0 * R * Λ_abs * R_inv * (u_rr - u_ll))
 
     return SVector(diss[1], diss[2], diss[3])
 end
@@ -506,11 +513,11 @@ end
 
     theta = rho_f * abs(shear_stress(u, equations)) / (gravity * (rho_s - rho_f) * d_s)  # Shields stress
 
-    Q = d_s * sqrt(gravity * (rho_s / rho_f - 1.0) * d_s) # Characteristic discharge
+    Q = d_s * sqrt(gravity * (rho_s / rho_f - 1) * d_s) # Characteristic discharge
 
     return (porosity_inv * Q * sign(theta) * k_1 * theta^m_1 *
-            (max(theta - k_2 * theta_c, 0.0))^m_2 *
-            (max(sqrt(theta) - k_3 * sqrt(theta_c), 0.0))^m_3)
+            (max(theta - k_2 * theta_c, 0))^m_2 *
+            (max(sqrt(theta) - k_3 * sqrt(theta_c), 0))^m_3)
 end
 
 # Compute the sediment discharge for the Grass model
@@ -553,7 +560,7 @@ end
 
     v = velocity(u, equations)
 
-    w1 = r * (gravity * (h + h_b) - 0.5 * v^2)
+    w1 = r * (gravity * (h + h_b) - 0.5f0 * v^2)
     w2 = r * v
     w3 = gravity * (r * h + h_b)
 
@@ -585,7 +592,8 @@ end
     v = velocity(u, equations)
     (; gravity, r) = equations
 
-    return 0.5 * r * h * v^2 + 0.5 * gravity * (r * h^2 + h_b^2) + r * gravity * h * h_b
+    return 0.5f0 * r * h * v^2 + 0.5f0 * gravity * (r * h^2 + h_b^2) +
+           r * gravity * h * h_b
 end
 
 # Calculate the error for the "lake-at-rest" test case where H = h + h_b should
@@ -605,16 +613,16 @@ end
     r = equations.r
 
     # Workaround to avoid division by zero, when computing the effective sediment height
-    if v > eps()
+    if abs(v) > eps(typeof(h))
         h_s = q_s(u, equations) / v
         # Compute gradients of q_s using automatic differentiation.
         # Introduces a closure to make q_s a function of u only. This is necessary since the
         # gradient function only accepts functions of one variable.
         dq_s_dh, dq_s_dhv, _ = Trixi.ForwardDiff.gradient(u -> q_s(u, equations), u)
     else
-        h_s = 0.0
-        dq_s_dh = 0.0
-        dq_s_dhv = 0.0
+        h_s = 0
+        dq_s_dh = 0
+        dq_s_dhv = 0
     end
 
     # Coefficient for the original cubic equation ax^3 + bx^2 + cx + d

test/test_tree_1d.jl

@@ -1074,14 +1074,14 @@ end # MLSWE
         @test_trixi_include(joinpath(EXAMPLES_DIR,
                                      "elixir_shallowwater_exner_channel.jl"),
                             l2=[
-                                0.061254947613784645,
-                                0.0042920880585939165,
-                                0.06170746499938789
+                                0.0612549484022613,
+                                0.004292092111911898,
+                                0.06170746566027052
                             ],
                             linf=[
-                                0.5552555774807875,
-                                0.009352028888004682,
-                                0.549962205546136
+                                0.555255659544283,
+                                0.009352017074210295,
+                                0.5499622869285822
                             ])
         # Ensure that we do not have excessive memory allocations
         # (e.g., from type instabilities)
@@ -1117,6 +1117,29 @@ end # MLSWE
             @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
         end
     end
+
+    @trixi_testset "elixir_shallowwater_exner_dam_break_symmetric.jl" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                     "elixir_shallowwater_exner_dam_break_symmetric.jl"),
+                            l2=[
+                                0.2965046413567244,
+                                0.41070730713108367,
+                                0.025179300052016823
+                            ],
+                            linf=[
+                                0.7785005894122902,
+                                0.9193314051729015,
+                                0.06639643395379602
+                            ])
+        # Ensure that we do not have excessive memory allocations
+        # (e.g., from type instabilities)
+        let
+            t = sol.t[end]
+            u_ode = sol.u[end]
+            du_ode = similar(u_ode)
+            @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        end
+    end
 end # SWE-Exner
 end # TreeMesh1D