rework thresholds

src/equations/equations.jl

@@ -41,4 +41,12 @@ Retrieve the number of layers from an equation instance of the `AbstractShallowW
                                                                                     }
     NLAYERS
 end
+
+# Provide default thresholds dependend on number format (Currently default thresholds are only provided
+# for Float64)
+default_threshold_partially_wet(::Type{Float64}) = 1e-4
+default_threshold_partially_wet(catchall) = throw(ArgumentError("threshold_partially_wet must be provided for non-Float64 types"))
+
+default_threshold_desingularization(::Type{Float64}) = 1e-10
+default_threshold_desingularization(catchall) = throw(ArgumentError("threshold_desingularization must be provided for non-Float64 types"))
 end # @muladd

src/equations/shallow_water_multilayer_1d.jl

@@ -30,11 +30,15 @@ nonconservative term, which has some benefits for the design of well-balanced sc
 The additional quantity ``H_0`` is also available to store a reference value for the total water
 height that is useful to set initial conditions or test the "lake-at-rest" well-balancedness.
 
-Also, there are two thresholds which prevent numerical problems as well as instabilities. Both of 
-them do not have to be passed, as default values are defined within the struct. The limiters are 
+Also, there are two thresholds which prevent numerical problems as well as instabilities. The limiters are 
 used in [`PositivityPreservingLimiterShallowWater`](@ref) on the water height. `threshold_limiter` 
 acts as a (small) shift on the initial condition and cutoff before the next time step, whereas 
-`threshold_desingularization` is used in the velocity desingularization.
+`threshold_desingularization` is used in the velocity desingularization. A third 
+`threshold_partially_wet` is applied on the water height to define "partially wet" elements in 
+[`IndicatorHennemannGassnerShallowWater`](@ref), that are then calculated with a pure FV method to
+ensure well-balancedness. For `Float64` no threshold needs to be passed, as default values are 
+defined within the struct. For other number formats `threshold_desingularization` and `threshold_partially_wet` 
+must be provided.
 
 The bottom topography function ``b(x)`` is set inside the initial condition routine
 for a particular problem setup.
@@ -62,12 +66,14 @@ struct ShallowWaterMultiLayerEquations1D{NVARS, NLAYERS, RealT <: Real} <:
     H0::RealT        # constant "lake-at-rest" total water height
     threshold_limiter::RealT    # threshold for the positivity-limiter
     threshold_desingularization::RealT  # threshold for velocity desingularization
+    threshold_partially_wet::RealT  # threshold to define partially wet elements
     rhos::SVector{NLAYERS, RealT} # Vector of layer densities
 
     function ShallowWaterMultiLayerEquations1D{NVARS, NLAYERS, RealT}(gravity::RealT,
                                                                       H0::RealT,
                                                                       threshold_limiter::RealT,
                                                                       threshold_desingularization::RealT,
+                                                                      threshold_partially_wet::RealT,
                                                                       rhos::SVector{NLAYERS,
                                                                                     RealT}) where {
                                                                                                    NVARS,
@@ -80,7 +86,8 @@ struct ShallowWaterMultiLayerEquations1D{NVARS, NLAYERS, RealT <: Real} <:
             throw(ArgumentError("densities must be in increasing order (rhos[1] < rhos[2] < ... < rhos[NLAYERS])"))
         min(rhos...) > 0 || throw(ArgumentError("densities must be positive"))
 
-        new(gravity, H0, threshold_limiter, threshold_desingularization, rhos)
+        new(gravity, H0, threshold_limiter, threshold_desingularization,
+            threshold_partially_wet, rhos)
     end
 end
 
@@ -91,7 +98,9 @@ end
 function ShallowWaterMultiLayerEquations1D(; gravity_constant,
                                            H0 = zero(gravity_constant),
                                            threshold_limiter = nothing,
-                                           threshold_desingularization = nothing, rhos)
+                                           threshold_desingularization = nothing,
+                                           threshold_partially_wet = nothing,
+                                           rhos)
 
     # Promote all variables to a common type
     _rhos = promote(rhos...)
@@ -103,7 +112,10 @@ function ShallowWaterMultiLayerEquations1D(; gravity_constant,
         threshold_limiter = 5 * eps(RealT)
     end
     if threshold_desingularization === nothing
-        threshold_desingularization = 1e-10
+        threshold_desingularization = default_threshold_desingularization(RealT)
+    end
+    if threshold_partially_wet === nothing
+        threshold_partially_wet = default_threshold_partially_wet(RealT)
     end
 
     # Extract number of layers and variables
@@ -114,6 +126,7 @@ function ShallowWaterMultiLayerEquations1D(; gravity_constant,
                                                                     H0,
                                                                     threshold_limiter,
                                                                     threshold_desingularization,
+                                                                    threshold_partially_wet,
                                                                     __rhos)
 end
 

src/equations/shallow_water_wet_dry_1d.jl

@@ -28,7 +28,11 @@ Also, there are two thresholds which prevent numerical problems as well as insta
 have to be passed, as default values are defined within the struct. The first one, `threshold_limiter`, is
 used in [`PositivityPreservingLimiterShallowWater`](@ref) on the water height, as a (small) shift on the initial
 condition and cutoff before the next time step. The second one, `threshold_wet`, is applied on the water height to
-define when the flow is "wet" before calculating the numerical flux.
+define when the flow is "wet" before calculating the numerical flux. A third 
+`threshold_partially_wet` is applied on the water height to define "partially wet" elements in 
+[`IndicatorHennemannGassnerShallowWater`](@ref), that are then calculated with a pure FV method to
+ensure well-balancedness. For `Float64` no threshold needs to be passed, as default values are 
+defined within the struct. For other number formats `threshold_partially_wet` must be provided.
 
 The bottom topography function ``b(x)`` is set inside the initial condition routine
 for a particular problem setup. To test the conservative form of the SWE one can set the bottom topography
@@ -62,6 +66,10 @@ struct ShallowWaterEquationsWetDry1D{RealT <: Real} <:
     # before calculating the numerical flux.
     # Default is 5*eps() which in double precision is ≈1e-15.
     threshold_wet::RealT
+    # `threshold_partially_wet` used in `IndicatorHennemannGassnerShallowWater` on the water height 
+    # to define "partially wet" elements. Those elements are calculated with a pure FV method to 
+    # ensure well-balancedness. Default in double precision is 1e-4. 
+    threshold_partially_wet::RealT
     # Standard shallow water equations for dispatch on Trixi.jl functions 
     basic_swe::ShallowWaterEquations1D{RealT}
 end
@@ -73,21 +81,25 @@ end
 # Strict default values for thresholds that performed well in many numerical experiments
 function ShallowWaterEquationsWetDry1D(; gravity_constant, H0 = zero(gravity_constant),
                                        threshold_limiter = nothing,
-                                       threshold_wet = nothing)
+                                       threshold_wet = nothing,
+                                       threshold_partially_wet = nothing)
     T = promote_type(typeof(gravity_constant), typeof(H0))
     if threshold_limiter === nothing
         threshold_limiter = 500 * eps(T)
     end
     if threshold_wet === nothing
         threshold_wet = 5 * eps(T)
     end
+    if threshold_partially_wet === nothing
+        threshold_partially_wet = default_threshold_partially_wet(T)
+    end
     # Construct the standard SWE for dispatch. Even though the `basic_swe` already store the 
     # gravity constant and the total water height, we store an extra copy in 
     # `ShallowWaterEquationsWetDry1D` for convenience.
     basic_swe = ShallowWaterEquations1D(gravity_constant = gravity_constant, H0 = H0)
 
     ShallowWaterEquationsWetDry1D(gravity_constant, H0, threshold_limiter,
-                                  threshold_wet, basic_swe)
+                                  threshold_wet, threshold_partially_wet, basic_swe)
 end
 
 Trixi.have_nonconservative_terms(::ShallowWaterEquationsWetDry1D) = True()

src/equations/shallow_water_wet_dry_2d.jl

@@ -31,7 +31,11 @@ Also, there are two thresholds which prevent numerical problems as well as insta
 have to be passed, as default values are defined within the struct. The first one, `threshold_limiter`, is
 used in [`PositivityPreservingLimiterShallowWater`](@ref) on the water height, as a (small) shift on the initial
 condition and cutoff before the next time step. The second one, `threshold_wet`, is applied on the water height to
-define when the flow is "wet" before calculating the numerical flux.
+define when the flow is "wet" before calculating the numerical flux. A third 
+`threshold_partially_wet` is applied on the water height to define "partially wet" elements in 
+[`IndicatorHennemannGassnerShallowWater`](@ref), that are then calculated with a pure FV method to
+ensure well-balancedness. For `Float64` no threshold needs to be passed, as default values are 
+defined within the struct. For other number formats  `threshold_partially_wet` must be provided.
 
 The bottom topography function ``b(x,y)`` is set inside the initial condition routine
 for a particular problem setup. To test the conservative form of the SWE one can set the bottom topography
@@ -65,6 +69,10 @@ struct ShallowWaterEquationsWetDry2D{RealT <: Real} <:
     # before calculating the numerical flux.
     # Default is 5*eps() which in double precision is ≈1e-15.
     threshold_wet::RealT
+    # `threshold_partially_wet` used in `IndicatorHennemannGassnerShallowWater` on the water height 
+    # to define "partially wet" elements. Those elements are calculated with a pure FV method to 
+    # ensure well-balancedness. Default in double precision is 1e-4. 
+    threshold_partially_wet::RealT
     # Standard shallow water equations for dispatch on Trixi.jl functions 
     basic_swe::ShallowWaterEquations2D{RealT}
 end
@@ -76,21 +84,25 @@ end
 # Strict default values for thresholds that performed well in many numerical experiments
 function ShallowWaterEquationsWetDry2D(; gravity_constant, H0 = zero(gravity_constant),
                                        threshold_limiter = nothing,
-                                       threshold_wet = nothing)
+                                       threshold_wet = nothing,
+                                       threshold_partially_wet = nothing)
     T = promote_type(typeof(gravity_constant), typeof(H0))
     if threshold_limiter === nothing
         threshold_limiter = 500 * eps(T)
     end
     if threshold_wet === nothing
         threshold_wet = 5 * eps(T)
     end
+    if threshold_partially_wet === nothing
+        threshold_partially_wet = default_threshold_partially_wet(T)
+    end
     # Construct the standard SWE for dispatch. Even though the `basic_swe` already store the 
     # gravity constant and the total water height, we store an extra copy in 
     # `ShallowWaterEquationsWetDry2D` for convenience.
     basic_swe = ShallowWaterEquations2D(gravity_constant = gravity_constant, H0 = H0)
 
     ShallowWaterEquationsWetDry2D(gravity_constant, H0, threshold_limiter,
-                                  threshold_wet, basic_swe)
+                                  threshold_wet, threshold_partially_wet, basic_swe)
 end
 
 Trixi.have_nonconservative_terms(::ShallowWaterEquationsWetDry2D) = True()

src/solvers/indicators_1d.jl

@@ -45,7 +45,7 @@ function (indicator_hg::IndicatorHennemannGassnerShallowWater)(u::AbstractArray{
     can only be proven for the FV method (see Chen and Noelle).
     Therefore we set alpha to one regardless of its given maximum value. 
     =#
-    threshold_partially_wet = 1e-4
+    threshold_partially_wet = equations.threshold_partially_wet
 
     Trixi.@threaded for element in eachelement(dg, cache)
         indicator = indicator_threaded[Threads.threadid()]

src/solvers/indicators_2d.jl

@@ -42,7 +42,7 @@ function (indicator_hg::IndicatorHennemannGassnerShallowWater)(u::AbstractArray{
     # Well-balancedness of the scheme on partially wet elements with hydrostatic reconstruction
     # can only be proven for the FV method (see Chen and Noelle).
     # Therefore we set alpha to be one regardless of its given value from the modal indicator.
-    threshold_partially_wet = 1e-4
+    threshold_partially_wet = equations.threshold_partially_wet
 
     Trixi.@threaded for element in eachelement(dg, cache)
         indicator = indicator_threaded[Threads.threadid()]