address changes from code review

examples/tree_1d_dgsem/elixir_shallowwater_multilayer_beach.jl

@@ -6,6 +6,7 @@ using TrixiShallowWater
 ###############################################################################
 # Semidiscretization of the shallow water equations
 
+# By passing only a single value for rhos, the system recovers the standard shallow water equations
 equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 9.812, rhos = 1.0)
 
 """

examples/tree_1d_dgsem/elixir_shallowwater_multilayer_parabolic_bowl.jl

@@ -6,6 +6,7 @@ using TrixiShallowWater
 ###############################################################################
 # Semidiscretization of the shallow water equations
 
+# By passing only a single value for rhos, the system recovers the standard shallow water equations
 equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 9.81, rhos = 1.0)
 
 """

examples/tree_1d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry.jl

@@ -75,8 +75,8 @@ function initial_condition_twolayer_well_balanced_wet_dry(perturbation::Bool,
         end
 
         # Set zero initial velocity
-        v1_upper = 0.0
-        v1_lower = 0.0
+        v1_upper = zero(H_upper)
+        v1_lower = zero(H_upper)
 
         #= 
         It is mandatory to shift the water level at dry areas to make sure the water height h
@@ -140,34 +140,6 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
 tspan = (0.0, 10.0)
 ode = semidiscretize(semi, tspan)
 
-#################################################################################
-# Workaround to set a discontinuous water and bottom topography for
-# debugging and testing. Essentially, this is a slight augmentation of the
-# `compute_coefficients` where the `x` node value passed here is slightly
-# perturbed to the left / right in order to set a true discontinuity that avoids
-# the doubled value of the LGL nodes at a particular element interface.
-#
-# Note! The errors from the analysis callback are not important but the error
-# for this lake at rest test case `∑|H0-(h+b)|` should be near machine roundoff.
-
-# point to the data we want to augment
-u = Trixi.wrap_array(ode.u0, semi)
-# reset the initial condition
-for element in eachelement(semi.solver, semi.cache)
-    for i in eachnode(semi.solver)
-        x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations,
-                                       semi.solver, i, element)
-        # We know that the discontinuity is a vertical line. Slightly augment the x value by a factor
-        # of unit roundoff to avoid the repeted value from the LGL nodes at at interface.
-        if i == 1
-            x_node = SVector(nextfloat(x_node[1]))
-        elseif i == nnodes(semi.solver)
-            x_node = SVector(prevfloat(x_node[1]))
-        end
-        u_node = initial_condition(x_node, first(tspan), equations)
-        Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, element)
-    end
-end
 #############################################################################################
 
 summary_callback = SummaryCallback()
@@ -193,7 +165,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, sav
 ###############################################################################
 # run the simulation
 
-# use a Runge-Kutta method with automatic (error based) time step size control
+# use a Runge-Kutta method with CFL-based time step
 sol = solve(ode, SSPRK43(stage_limiter!);
             ode_default_options()..., callback = callbacks, adaptive = false, dt = 1.0);
 summary_callback() # print the timer summary

src/callbacks_stage/positivity_shallow_water.jl

@@ -8,9 +8,6 @@
 """
     PositivityPreservingLimiterShallowWater(; variables)
 
-!!! warning "Experimental code"
-    This is an experimental feature and may change in future releases.
-
 The limiter is specifically designed for the shallow water equations.
 It is applied to all scalar `variables` in their given order
 using the defined `threshold_limiter` from the equations struct 
@@ -50,6 +47,10 @@ The specific implementation for the [`ShallowWaterMultiLayerEquations1D](@ref) i
   Positivity-preserving well-balanced discontinuous Galerkin methods for the shallow water equations
   on unstructured triangular meshes
   [doi: 10.1007/s10915-012-9644-4](https://doi.org/10.1007/s10915-012-9644-4)
+
+
+!!! warning "Experimental code"
+    This is an experimental feature and may change in future releases.
 """
 struct PositivityPreservingLimiterShallowWater{N, Variables <: NTuple{N, Any}}
     variables::Variables

src/callbacks_stage/positivity_shallow_water_dg1d.jl

@@ -87,12 +87,12 @@ function limiter_shallow_water!(u, threshold::Real, variable,
     return nothing
 end
 
-# !!! warning "Experimental code"
-#     This is an experimental feature and may change in future releases.
-
 # Note that for the `ShallowWaterMultiLayerEquations1D` only the waterheight `h` is limited in
 # each layer. Furthermore, a velocity desingularization is applied after the limiting to avoid
 # numerical problems near dry states.
+#
+# !!! warning "Experimental code"
+#     This is an experimental feature and may change in future releases.
 function limiter_shallow_water!(u, threshold::Real, variable,
                                 mesh::Trixi.AbstractMesh{1},
                                 equations::ShallowWaterMultiLayerEquations1D,

src/equations/shallow_water_multilayer_1d.jl

@@ -135,13 +135,13 @@ function Trixi.varnames(::typeof(cons2cons),
     return (waterheight..., momentum..., "b")
 end
 
-# We use the interface heights, H = ∑h + b as primitive variables for easier visualization and setting initial
+# We use the total layer heights, H = ∑h + b as primitive variables for easier visualization and setting initial
 # conditions
 function Trixi.varnames(::typeof(cons2prim),
                         equations::ShallowWaterMultiLayerEquations1D)
-    interface_height = ntuple(n -> "H" * string(n), Val(nlayers(equations)))
+    total_layer_height = ntuple(n -> "H" * string(n), Val(nlayers(equations)))
     velocity = ntuple(n -> "v" * string(n) * "_1", Val(nlayers(equations)))
-    return (interface_height..., velocity..., "b")
+    return (total_layer_height..., velocity..., "b")
 end
 
 # Set initial conditions at physical location `x` for time `t`
@@ -379,9 +379,6 @@ end
 """
     hydrostatic_reconstruction_ersing_etal(u_ll, u_rr, equations::ShallowWaterMultiLayerEquations1D)
 
-!!! warning "Experimental code"
-    This is an experimental feature and may change in future releases.
-
 A particular type of hydrostatic reconstruction of the water height and bottom topography to 
 guarantee well-balancedness in the presence of wet/dry transitions and entropy stability for the 
 [`ShallowWaterMultiLayerEquations1D`](@ref). 
@@ -390,6 +387,9 @@ surface numerical flux at the interface. The key idea is a piecewise linear reco
 bottom topography and water height interfaces using subcells, where the bottom topography is allowed 
 to be discontinuous. 
 Use in combination with the generic numerical flux routine [`Trixi.FluxHydrostaticReconstruction`](@extref).
+
+!!! warning "Experimental code"
+    This is an experimental feature and may change in future releases.
 """
 @inline function hydrostatic_reconstruction_ersing_etal(u_ll, u_rr,
                                                         equations::ShallowWaterMultiLayerEquations1D)
@@ -415,15 +415,15 @@ Use in combination with the generic numerical flux routine [`Trixi.FluxHydrostat
         end
     end
 
-    # Calculate interface heights
+    # Calculate total layer heights
     H_ll = waterheight(cons2prim(u_ll, equations), equations)
     H_rr = waterheight(cons2prim(u_rr, equations), equations)
 
     # Discontinuous reconstruction of the bottom topography
     b_ll_star = min(H_ll[1], max(b_rr, b_ll))
     b_rr_star = min(H_rr[1], max(b_rr, b_ll))
 
-    # Calculate reconstructed interface heights
+    # Calculate reconstructed total layer heights
     H_ll_star = max.(H_ll, b_ll_star)
     H_rr_star = max.(H_rr, b_rr_star)
 
@@ -501,7 +501,7 @@ end
     h = waterheight(u, equations)
     b = u[end]
 
-    # Initialize interface height
+    # Initialize total layer height
     H = MVector{nlayers(equations), real(equations)}(undef)
     for i in reverse(eachlayer(equations))
         if i == nlayers(equations)
@@ -517,9 +517,9 @@ end
 
 # Convert primitive to conservative variables
 @inline function Trixi.prim2cons(prim, equations::ShallowWaterMultiLayerEquations1D)
-    # To extract the interface height and velocity we reuse the waterheight and momentum functions 
+    # To extract the total layer height and velocity we reuse the waterheight and momentum functions 
     # from the conservative variables.
-    H = waterheight(prim, equations)    # For primitive variables this extracts the interface height
+    H = waterheight(prim, equations)    # For primitive variables this extracts the total layer height
     v = momentum(prim, equations)       # For primitive variables this extracts the velocity
     b = prim[end]
 
@@ -556,9 +556,9 @@ end
 
     # Calculate entropy variables in each layer
     for i in eachlayer(equations)
-        # Compute w1[i] = ρ[i] * g * (b + ∑h[k] + ∑σ[k] * h[k]), where σ[k] = ρ[k] / ρ[i] denotes 
-        # the density ratio of different layers
-        w1 = equations.rhos[i] * (g * b - 0.5 * v[i]^2)
+        # Compute w1[i] = ρ[i] * g * (b + ∑h[k] + ∑σ[k] * h[k]) - 0.5 * ρ[i] * v[i]^2, 
+        # where σ[k] = ρ[k] / ρ[i] denotes the density ratio of different layers
+        w1 = equations.rhos[i] * (g * b) - 0.5 * equations.rhos[i] * v[i]^2
         for j in eachlayer(equations)
             if j < i
                 w1 += equations.rhos[i] * g *