From 3adf848ac4d4d70f1f1e87143d0c719101c14ab3 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 11 Dec 2024 12:26:44 +0100 Subject: [PATCH] Thread Parallel Reduction for `integrate_via_indices` --- src/callbacks_step/analysis_dg1d.jl | 27 ++++++++++---------- src/callbacks_step/analysis_dg2d.jl | 5 ++-- src/callbacks_step/analysis_dg2d_parallel.jl | 2 +- src/callbacks_step/analysis_dg3d.jl | 5 ++-- src/callbacks_step/analysis_dg3d_parallel.jl | 2 +- src/solvers/fdsbp_tree/fdsbp_1d.jl | 2 +- src/solvers/fdsbp_tree/fdsbp_2d.jl | 2 +- src/solvers/fdsbp_tree/fdsbp_3d.jl | 2 +- src/solvers/fdsbp_unstructured/fdsbp_2d.jl | 3 ++- 9 files changed, 27 insertions(+), 23 deletions(-) diff --git a/src/callbacks_step/analysis_dg1d.jl b/src/callbacks_step/analysis_dg1d.jl index b4d49a39d5a..82064c7cb85 100644 --- a/src/callbacks_step/analysis_dg1d.jl +++ b/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 diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 628c3071cbc..d8da0ec7f70 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/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] * diff --git a/src/callbacks_step/analysis_dg2d_parallel.jl b/src/callbacks_step/analysis_dg2d_parallel.jl index 10b6e30f95c..ba2989f721d 100644 --- a/src/callbacks_step/analysis_dg2d_parallel.jl +++ b/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] * diff --git a/src/callbacks_step/analysis_dg3d.jl b/src/callbacks_step/analysis_dg3d.jl index dfd5696b5fb..8bb6c630162 100644 --- a/src/callbacks_step/analysis_dg3d.jl +++ b/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] * diff --git a/src/callbacks_step/analysis_dg3d_parallel.jl b/src/callbacks_step/analysis_dg3d_parallel.jl index ab9ba6a0055..23aaa7319bc 100644 --- a/src/callbacks_step/analysis_dg3d_parallel.jl +++ b/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] * diff --git a/src/solvers/fdsbp_tree/fdsbp_1d.jl b/src/solvers/fdsbp_tree/fdsbp_1d.jl index 0de0cff4851..9017da591cc 100644 --- a/src/solvers/fdsbp_tree/fdsbp_1d.jl +++ b/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] * diff --git a/src/solvers/fdsbp_tree/fdsbp_2d.jl b/src/solvers/fdsbp_tree/fdsbp_2d.jl index 36afbbc022f..0d904a7ef39 100644 --- a/src/solvers/fdsbp_tree/fdsbp_2d.jl +++ b/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] * diff --git a/src/solvers/fdsbp_tree/fdsbp_3d.jl b/src/solvers/fdsbp_tree/fdsbp_3d.jl index 0c3f18b6d6e..e212423701f 100644 --- a/src/solvers/fdsbp_tree/fdsbp_3d.jl +++ b/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] * diff --git a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl index cbe11ac6ac9..7e66e0d472b 100644 --- a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl +++ b/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] *