Mean value in Zhang-Shu limiter on curved meshes (#1945)

* solution mean value computation must take metric terms into account on mesh types other than TreeMesh

* introduce helper function to compute mean value on 2d element

* apply formatter

* introduce helper function in 3d

* simplify mean value computation. no need for a separate function

* fix setting total_volume to be zero

* put calc_mean_value functions back in

* apply formatter

* add generic get_inverse_jacobian function and remove other specializations

* swap order of element and indices for consistency

* revert ordering. see if tests actually run

* fix index order

* fix mistake in TreeMesh version

* clarify the tree mesh version

src/callbacks_stage/positivity_zhang_shu_dg2d.jl

@@ -8,6 +8,7 @@
 function limiter_zhang_shu!(u, threshold::Real, variable,
                             mesh::AbstractMesh{2}, equations, dg::DGSEM, cache)
     @unpack weights = dg.basis
+    @unpack inverse_jacobian = cache.elements
 
     @threaded for element in eachelement(dg, cache)
         # determine minimum value
@@ -22,12 +23,16 @@ function limiter_zhang_shu!(u, threshold::Real, variable,
 
         # compute mean value
         u_mean = zero(get_node_vars(u, equations, dg, 1, 1, element))
+        total_volume = zero(eltype(u))
         for j in eachnode(dg), i in eachnode(dg)
+            volume_jacobian = abs(inv(get_inverse_jacobian(inverse_jacobian, mesh,
+                                                           i, j, element)))
             u_node = get_node_vars(u, equations, dg, i, j, element)
-            u_mean += u_node * weights[i] * weights[j]
+            u_mean += u_node * weights[i] * weights[j] * volume_jacobian
+            total_volume += weights[i] * weights[j] * volume_jacobian
         end
-        # note that the reference element is [-1,1]^ndims(dg), thus the weights sum to 2
-        u_mean = u_mean / 2^ndims(mesh)
+        # 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).

src/callbacks_stage/positivity_zhang_shu_dg3d.jl

@@ -8,6 +8,7 @@
 function limiter_zhang_shu!(u, threshold::Real, variable,
                             mesh::AbstractMesh{3}, equations, dg::DGSEM, cache)
     @unpack weights = dg.basis
+    @unpack inverse_jacobian = cache.elements
 
     @threaded for element in eachelement(dg, cache)
         # determine minimum value
@@ -22,12 +23,16 @@ function limiter_zhang_shu!(u, threshold::Real, variable,
 
         # compute mean value
         u_mean = zero(get_node_vars(u, equations, dg, 1, 1, 1, element))
+        total_volume = zero(eltype(u))
         for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+            volume_jacobian = abs(inv(get_inverse_jacobian(inverse_jacobian, mesh,
+                                                           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]
+            u_mean += u_node * weights[i] * weights[j] * weights[k] * volume_jacobian
+            total_volume += weights[i] * weights[j] * weights[k] * volume_jacobian
         end
-        # note that the reference element is [-1,1]^ndims(dg), thus the weights sum to 2
-        u_mean = u_mean / 2^ndims(mesh)
+        # 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).

src/solvers/dgsem_structured/dg.jl

@@ -76,6 +76,14 @@ end
     end
 end
 
+@inline function get_inverse_jacobian(inverse_jacobian,
+                                      mesh::Union{StructuredMesh, StructuredMeshView,
+                                                  UnstructuredMesh2D, P4estMesh,
+                                                  T8codeMesh},
+                                      indices...)
+    return inverse_jacobian[indices...]
+end
+
 include("containers.jl")
 include("dg_1d.jl")
 include("dg_2d.jl")

src/solvers/dgsem_tree/dg.jl

@@ -42,6 +42,12 @@ function volume_jacobian(element, mesh::TreeMesh, cache)
     return inv(cache.elements.inverse_jacobian[element])^ndims(mesh)
 end
 
+@inline function get_inverse_jacobian(inverse_jacobian, mesh::TreeMesh,
+                                      indices...)
+    element = last(indices)
+    return inverse_jacobian[element]
+end
+
 # Indicators used for shock-capturing and AMR
 include("indicators.jl")
 include("indicators_1d.jl")