Mean value in Zhang-Shu limiter on curved meshes

src/callbacks_stage/positivity_zhang_shu_dg2d.jl

@@ -6,7 +6,7 @@
 #! format: noindent
 
 function limiter_zhang_shu!(u, threshold::Real, variable,
-                            mesh::AbstractMesh{2}, equations, dg::DGSEM, cache)
+                            mesh::TreeMesh{2}, equations, dg::DGSEM, cache)
     @unpack weights = dg.basis
 
     @threaded for element in eachelement(dg, cache)
@@ -42,4 +42,48 @@ function limiter_zhang_shu!(u, threshold::Real, variable,
 
     return nothing
 end
+
+function limiter_zhang_shu!(u, threshold::Real, variable,
+                            mesh::Union{StructuredMesh{2}, StructuredMeshView{2},
+                                        UnstructuredMesh2D, P4estMesh{2},
+                                        T8codeMesh{2}},
+                            equations, dg::DGSEM, cache)
+    @unpack weights = dg.basis
+
+    @threaded for element in eachelement(dg, cache)
+        # determine minimum value
+        value_min = typemax(eltype(u))
+        for j in eachnode(dg), i in eachnode(dg)
+            u_node = get_node_vars(u, equations, dg, i, j, element)
+            value_min = min(value_min, variable(u_node, equations))
+        end
+
+        # detect if limiting is necessary
+        value_min < threshold || continue
+
+        # compute mean value
+        u_mean = zero(get_node_vars(u, equations, dg, 1, 1, element))
+        total_volume = zero(real(mesh))
+        for j in eachnode(dg), i in eachnode(dg)
+            jacobian_node = abs(inv(cache.elements.inverse_jacobian[i, j, element]))
+            u_node = get_node_vars(u, equations, dg, i, j, element)
+            u_mean += u_node * weights[i] * weights[j] * jacobian_node
+            total_volume += weights[i] * weights[j] * jacobian_node
+        end
+        # normalize with the total volume
+        u_mean = u_mean / total_volume
+
+        # We compute the value directly with the mean values, as we assume that
+        # Jensen's inequality holds (e.g. pressure for compressible Euler equations).
+        value_mean = variable(u_mean, equations)
+        theta = (value_mean - threshold) / (value_mean - value_min)
+        for j in eachnode(dg), i in eachnode(dg)
+            u_node = get_node_vars(u, equations, dg, i, j, element)
+            set_node_vars!(u, theta * u_node + (1 - theta) * u_mean,
+                           equations, dg, i, j, element)
+        end
+    end
+
+    return nothing
+end
 end # @muladd

src/callbacks_stage/positivity_zhang_shu_dg3d.jl

@@ -6,7 +6,7 @@
 #! format: noindent
 
 function limiter_zhang_shu!(u, threshold::Real, variable,
-                            mesh::AbstractMesh{3}, equations, dg::DGSEM, cache)
+                            mesh::TreeMesh{3}, equations, dg::DGSEM, cache)
     @unpack weights = dg.basis
 
     @threaded for element in eachelement(dg, cache)
@@ -42,4 +42,46 @@ function limiter_zhang_shu!(u, threshold::Real, variable,
 
     return nothing
 end
+
+function limiter_zhang_shu!(u, threshold::Real, variable,
+                            mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}},
+                            equations, dg::DGSEM, cache)
+    @unpack weights = dg.basis
+
+    @threaded for element in eachelement(dg, cache)
+        # determine minimum value
+        value_min = typemax(eltype(u))
+        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+            u_node = get_node_vars(u, equations, dg, i, j, k, element)
+            value_min = min(value_min, variable(u_node, equations))
+        end
+
+        # detect if limiting is necessary
+        value_min < threshold || continue
+
+        # compute mean value
+        u_mean = zero(get_node_vars(u, equations, dg, 1, 1, 1, element))
+        total_volume = zero(real(mesh))
+        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+            jacobian_node = abs(inv(cache.elements.inverse_jacobian[i, j, k, element]))
+            u_node = get_node_vars(u, equations, dg, i, j, k, element)
+            u_mean += u_node * weights[i] * weights[j] * weights[k] * jacobian_node
+            total_volume += weights[i] * weights[j] * weights[k] * jacobian_node
+        end
+        # normalize with the total volume
+        u_mean = u_mean / total_volume
+
+        # We compute the value directly with the mean values, as we assume that
+        # Jensen's inequality holds (e.g. pressure for compressible Euler equations).
+        value_mean = variable(u_mean, equations)
+        theta = (value_mean - threshold) / (value_mean - value_min)
+        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+            u_node = get_node_vars(u, equations, dg, i, j, k, element)
+            set_node_vars!(u, theta * u_node + (1 - theta) * u_mean,
+                           equations, dg, i, j, k, element)
+        end
+    end
+
+    return nothing
+end
 end # @muladd