Thread Parallel Reduction for integrate_via_indices

src/callbacks_step/analysis_dg1d.jl

@@ -120,51 +120,52 @@ function calc_error_norms(func, u, t, analyzer,
 end
 
 function integrate_via_indices(func::Func, u,
-                               mesh::StructuredMesh{1}, equations, dg::DGSEM, cache,
+                               mesh::TreeMesh{1}, equations, dg::DGSEM, cache,
                                args...; normalize = true) where {Func}
     @unpack weights = dg.basis
 
     # Initialize integral with zeros of the right shape
     integral = zero(func(u, 1, 1, equations, dg, args...))
-    total_volume = zero(real(mesh))
 
     # Use quadrature to numerically integrate over entire domain
-    for element in eachelement(dg, cache)
+    @batch reduction=(+, integral) for element in eachelement(dg, cache)
+        volume_jacobian_ = volume_jacobian(element, mesh, cache)
         for i in eachnode(dg)
-            jacobian_volume = abs(inv(cache.elements.inverse_jacobian[i, element]))
-            integral += jacobian_volume * weights[i] *
+            integral += volume_jacobian_ * weights[i] *
                         func(u, i, element, equations, dg, args...)
-            total_volume += jacobian_volume * weights[i]
         end
     end
+
     # Normalize with total volume
     if normalize
-        integral = integral / total_volume
+        integral = integral / total_volume(mesh)
     end
 
     return integral
 end
 
 function integrate_via_indices(func::Func, u,
-                               mesh::TreeMesh{1}, equations, dg::DGSEM, cache,
+                               mesh::StructuredMesh{1}, equations, dg::DGSEM, cache,
                                args...; normalize = true) where {Func}
     @unpack weights = dg.basis
 
     # Initialize integral with zeros of the right shape
     integral = zero(func(u, 1, 1, equations, dg, args...))
+    total_volume = zero(real(mesh))
 
     # Use quadrature to numerically integrate over entire domain
-    for element in eachelement(dg, cache)
-        volume_jacobian_ = volume_jacobian(element, mesh, cache)
+    @batch reduction=((+, integral), (+, total_volume)) for element in eachelement(dg,
+                                                                                   cache)
         for i in eachnode(dg)
-            integral += volume_jacobian_ * weights[i] *
+            jacobian_volume = abs(inv(cache.elements.inverse_jacobian[i, element]))
+            integral += jacobian_volume * weights[i] *
                         func(u, i, element, equations, dg, args...)
+            total_volume += jacobian_volume * weights[i]
         end
     end
-
     # Normalize with total volume
     if normalize
-        integral = integral / total_volume(mesh)
+        integral = integral / total_volume
     end
 
     return integral

src/callbacks_step/analysis_dg2d.jl

@@ -190,7 +190,7 @@ function integrate_via_indices(func::Func, u,
     integral = zero(func(u, 1, 1, 1, equations, dg, args...))
 
     # Use quadrature to numerically integrate over entire domain
-    for element in eachelement(dg, cache)
+    @batch reduction=(+, integral) for element in eachelement(dg, cache)
         volume_jacobian_ = volume_jacobian(element, mesh, cache)
         for j in eachnode(dg), i in eachnode(dg)
             integral += volume_jacobian_ * weights[i] * weights[j] *
@@ -219,7 +219,8 @@ function integrate_via_indices(func::Func, u,
     total_volume = zero(real(mesh))
 
     # Use quadrature to numerically integrate over entire domain
-    for element in eachelement(dg, cache)
+    @batch reduction=((+, integral), (+, total_volume)) for element in eachelement(dg,
+                                                                                   cache)
         for j in eachnode(dg), i in eachnode(dg)
             volume_jacobian = abs(inv(cache.elements.inverse_jacobian[i, j, element]))
             integral += volume_jacobian * weights[i] * weights[j] *

src/callbacks_step/analysis_dg2d_parallel.jl

@@ -187,7 +187,7 @@ function integrate_via_indices(func::Func, u,
     volume = zero(real(mesh))
 
     # Use quadrature to numerically integrate over entire domain
-    for element in eachelement(dg, cache)
+    @batch reduction=((+, integral), (+, volume)) for element in eachelement(dg, cache)
         for j in eachnode(dg), i in eachnode(dg)
             volume_jacobian = abs(inv(cache.elements.inverse_jacobian[i, j, element]))
             integral += volume_jacobian * weights[i] * weights[j] *

src/callbacks_step/analysis_dg3d.jl

@@ -218,7 +218,7 @@ function integrate_via_indices(func::Func, u,
     integral = zero(func(u, 1, 1, 1, 1, equations, dg, args...))
 
     # Use quadrature to numerically integrate over entire domain
-    for element in eachelement(dg, cache)
+    @batch reduction=(+, integral) for element in eachelement(dg, cache)
         volume_jacobian_ = volume_jacobian(element, mesh, cache)
         for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
             integral += volume_jacobian_ * weights[i] * weights[j] * weights[k] *
@@ -246,7 +246,8 @@ function integrate_via_indices(func::Func, u,
     total_volume = zero(real(mesh))
 
     # Use quadrature to numerically integrate over entire domain
-    for element in eachelement(dg, cache)
+    @batch reduction=((+, integral), (+, total_volume)) for element in eachelement(dg,
+                                                                                   cache)
         for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
             volume_jacobian = abs(inv(cache.elements.inverse_jacobian[i, j, k, element]))
             integral += volume_jacobian * weights[i] * weights[j] * weights[k] *

src/callbacks_step/analysis_dg3d_parallel.jl

@@ -82,7 +82,7 @@ function integrate_via_indices(func::Func, u,
     volume = zero(real(mesh))
 
     # Use quadrature to numerically integrate over entire domain
-    for element in eachelement(dg, cache)
+    @batch reduction=((+, integral), (+, volume)) for element in eachelement(dg, cache)
         for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
             volume_jacobian = abs(inv(cache.elements.inverse_jacobian[i, j, k, element]))
             integral += volume_jacobian * weights[i] * weights[j] * weights[k] *

src/solvers/fdsbp_tree/fdsbp_1d.jl

@@ -271,7 +271,7 @@ function integrate_via_indices(func::Func, u,
     integral = zero(func(u, 1, 1, equations, dg, args...))
 
     # Use quadrature to numerically integrate over entire domain
-    for element in eachelement(dg, cache)
+    @batch reduction=(+, integral) for element in eachelement(dg, cache)
         volume_jacobian_ = volume_jacobian(element, mesh, cache)
         for i in eachnode(dg)
             integral += volume_jacobian_ * weights[i] *

src/solvers/fdsbp_tree/fdsbp_2d.jl

@@ -326,7 +326,7 @@ function integrate_via_indices(func::Func, u,
     integral = zero(func(u, 1, 1, 1, equations, dg, args...))
 
     # Use quadrature to numerically integrate over entire domain
-    for element in eachelement(dg, cache)
+    @batch reduction=(+, integral) for element in eachelement(dg, cache)
         volume_jacobian_ = volume_jacobian(element, mesh, cache)
         for j in eachnode(dg), i in eachnode(dg)
             integral += volume_jacobian_ * weights[i] * weights[j] *

src/solvers/fdsbp_tree/fdsbp_3d.jl

@@ -378,7 +378,7 @@ function integrate_via_indices(func::Func, u,
     integral = zero(func(u, 1, 1, 1, 1, equations, dg, args...))
 
     # Use quadrature to numerically integrate over entire domain
-    for element in eachelement(dg, cache)
+    @batch reduction=(+, integral) for element in eachelement(dg, cache)
         volume_jacobian_ = volume_jacobian(element, mesh, cache)
         for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
             integral += volume_jacobian_ * weights[i] * weights[j] * weights[k] *

src/solvers/fdsbp_unstructured/fdsbp_2d.jl

@@ -247,7 +247,8 @@ function integrate_via_indices(func::Func, u,
     total_volume = zero(real(mesh))
 
     # Use quadrature to numerically integrate over entire domain
-    for element in eachelement(dg, cache)
+    @batch reduction=((+, integral), (+, total_volume)) for element in eachelement(dg,
+                                                                                   cache)
         for j in eachnode(dg), i in eachnode(dg)
             volume_jacobian = abs(inv(cache.elements.inverse_jacobian[i, j, element]))
             integral += volume_jacobian * weights[i] * weights[j] *