From 117489ee6423469a2dabe2d0ad95b69a955aed27 Mon Sep 17 00:00:00 2001
From: Patrick Ersing <114223904+patrickersing@users.noreply.github.com>
Date: Wed, 22 May 2024 08:40:24 +0200
Subject: [PATCH] add limiter fix for curvilinear meshes (#46)

---
 .../positivity_shallow_water_dg2d.jl          | 23 ++++++++++++++-----
 1 file changed, 17 insertions(+), 6 deletions(-)

diff --git a/src/callbacks_stage/positivity_shallow_water_dg2d.jl b/src/callbacks_stage/positivity_shallow_water_dg2d.jl
index 9acdc4d..bea301f 100644
--- a/src/callbacks_stage/positivity_shallow_water_dg2d.jl
+++ b/src/callbacks_stage/positivity_shallow_water_dg2d.jl
@@ -12,6 +12,7 @@ function limiter_shallow_water!(u, threshold::Real, variable,
                                 equations::ShallowWaterEquationsWetDry2D, dg::DGSEM,
                                 cache)
     @unpack weights = dg.basis
+    @unpack inverse_jacobian = cache.elements
 
     Trixi.@threaded for element in eachelement(dg, cache)
         # determine minimum value
@@ -26,12 +27,16 @@ function limiter_shallow_water!(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(Trixi.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).
@@ -96,6 +101,7 @@ function limiter_shallow_water!(u, threshold::Real, variable,
                                 equations::ShallowWaterMultiLayerEquations2D, dg::DGSEM,
                                 cache)
     @unpack weights = dg.basis
+    @unpack inverse_jacobian = cache.elements
 
     Trixi.@threaded for element in eachelement(dg, cache)
         # Limit layerwise
@@ -112,12 +118,17 @@ function limiter_shallow_water!(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(Trixi.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.
             # The waterheight `h` is limited independently in each layer.