From f28ca44ac9b22b26d7ef1ca9fc589df96b814724 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Tue, 5 Nov 2024 09:17:49 +0100 Subject: [PATCH] Add Simpson quadrature with trapzoidal extension Add an additional Simpson quadrature which uses trapzoid quadrature over one cell when the number of grid points is not compatible with the Simpson quadrature. Use this quadrature for the VparMu collision test. See merge request gysela-developpers/gyselalibxx!760 -------------------------------------------- --- .../collisions/test_collSpVparMu.cpp | 5 +- src/collisions/CollisionSpVparMu.hpp | 2 +- src/multipatch/connectivity/edge.hpp | 5 +- src/quadrature/simpson_quadrature.hpp | 94 +++++++++++++++++-- src/utils/geometry_descriptors.hpp | 4 + 5 files changed, 94 insertions(+), 16 deletions(-) create mode 100644 src/utils/geometry_descriptors.hpp diff --git a/simulations/geometryVparMu/collisions/test_collSpVparMu.cpp b/simulations/geometryVparMu/collisions/test_collSpVparMu.cpp index acd7ba620..b7d10eafd 100644 --- a/simulations/geometryVparMu/collisions/test_collSpVparMu.cpp +++ b/simulations/geometryVparMu/collisions/test_collSpVparMu.cpp @@ -22,6 +22,7 @@ #include "paraconfpp.hpp" #include "params.yaml.hpp" #include "pdi_out.yml.hpp" +#include "simpson_quadrature.hpp" #include "species_info.hpp" #include "species_init.hpp" @@ -90,8 +91,8 @@ int main(int argc, char** argv) DFieldMemVpar const coeff_intdvpar(neumann_spline_quadrature_coefficients< Kokkos::DefaultExecutionSpace>(idxrange_vpar, builder_vpar)); SplineMuBuilder const builder_mu(idxrange_mu); - DFieldMemMu const coeff_intdmu(neumann_spline_quadrature_coefficients< - Kokkos::DefaultExecutionSpace>(idxrange_mu, builder_mu)); + DFieldMemMu const coeff_intdmu(simpson_trapezoid_quadrature_coefficients_1d< + Kokkos::DefaultExecutionSpace>(idxrange_mu, Extremity::BACK)); CollisionInfo const collision_info(conf_collision); CollisionSpVparMu collision_operator( collision_info, diff --git a/src/collisions/CollisionSpVparMu.hpp b/src/collisions/CollisionSpVparMu.hpp index baa36df52..1938f1df0 100644 --- a/src/collisions/CollisionSpVparMu.hpp +++ b/src/collisions/CollisionSpVparMu.hpp @@ -172,11 +172,11 @@ class CollisionSpVparMu /* : public IRightHandSide */ : m_operator_handle {} , m_hat_As {"m_hat_As", collisions_dimensions::get_idx_range(fdistrib_idx_range)} , m_hat_Zs {"m_hat_Zs", collisions_dimensions::get_idx_range(fdistrib_idx_range)} + , m_coeff_AD {"m_coeff_AD", collisions_dimensions::get_idx_range(fdistrib_idx_range)} , m_mask_buffer_r {"m_mask_buffer_r", collisions_dimensions::get_idx_range(fdistrib_idx_range)} , m_mask_LIM {"m_mask_LIM", IdxRangeThetaR {collisions_dimensions::get_idx_range(fdistrib_idx_range)}} , m_B_norm {"m_B_norm", IdxRangeThetaR {collisions_dimensions::get_idx_range(fdistrib_idx_range)}} , m_Bstar_s {"m_Bstar_s", IdxRangeSpThetaRVpar {collisions_dimensions::get_idx_range(fdistrib_idx_range)}} - , m_coeff_AD {"m_coeff_AD", collisions_dimensions::get_idx_range(fdistrib_idx_range)} , m_mug {"m_mug", ddc::select(fdistrib_idx_range)} , m_vparg {"m_vparg", ddc::select(fdistrib_idx_range)} { diff --git a/src/multipatch/connectivity/edge.hpp b/src/multipatch/connectivity/edge.hpp index 40317c7dc..ab3eac8c7 100644 --- a/src/multipatch/connectivity/edge.hpp +++ b/src/multipatch/connectivity/edge.hpp @@ -1,12 +1,9 @@ // SPDX-License-Identifier: MIT - #pragma once #include #include "ddc_aliases.hpp" - - -enum Extremity { FRONT, BACK }; +#include "geometry_descriptors.hpp" /** diff --git a/src/quadrature/simpson_quadrature.hpp b/src/quadrature/simpson_quadrature.hpp index eb88bdca0..569d4e4e7 100644 --- a/src/quadrature/simpson_quadrature.hpp +++ b/src/quadrature/simpson_quadrature.hpp @@ -4,7 +4,9 @@ #include #include "ddc_aliases.hpp" +#include "geometry_descriptors.hpp" #include "quadrature_coeffs_nd.hpp" +#include "trapezoid_quadrature.hpp" /** * @brief Get the Simpson coefficients in 1D. @@ -14,15 +16,14 @@ * * (x_3-x_1)(2-(x_3-x_2)/(x_2-x_1)) / 6 f(x_1) + (x_3-x_1)^3 / 6 / (x_3-x_2) / (x_2-x_1) f(x_2) + (x_3-x_1)(2-(x_2-x_1)/(x_3-x_2)) / 6 f(x_3) * - * @param[in] idx_range - * The index range on which the quadrature will be carried out. - * - * @return The quadrature coefficients for the Simpson method defined on the provided index range. + * @param[in] coefficients + * The field where the quadrature coefficients should be saved. */ template -DFieldMem, typename ExecSpace::memory_space> simpson_quadrature_coefficients_1d( - IdxRange const& idx_range) +void fill_simpson_quadrature_coefficients_1d( + DField, typename ExecSpace::memory_space> coefficients) { + IdxRange idx_range = get_idx_range(coefficients); if constexpr (Grid1D::continuous_dimension_type::PERIODIC) { if (idx_range.size() % 2 != 0) { throw std::runtime_error( @@ -34,9 +35,6 @@ DFieldMem, typename ExecSpace::memory_space> simpson_quadrature "non-periodic direction."); } } - DFieldMem, typename ExecSpace::memory_space> coefficients_alloc(idx_range); - DField, typename ExecSpace::memory_space> const coefficients - = get_field(coefficients_alloc); Kokkos::parallel_for( "centre_left", @@ -92,10 +90,88 @@ DFieldMem, typename ExecSpace::memory_space> simpson_quadrature coefficients(idx_range.front()) += dx_sum * (2. * dx_r - dx_l) / (6. * dx_r); }); } +} +/** + * @brief Get the Simpson coefficients in 1D. + * + * Calculate the quadrature coefficients for the Simpson method defined on the provided index range. + * The non-uniform form of the Simpson quadrature is: + * + * (x_3-x_1)(2-(x_3-x_2)/(x_2-x_1)) / 6 f(x_1) + (x_3-x_1)^3 / 6 / (x_3-x_2) / (x_2-x_1) f(x_2) + (x_3-x_1)(2-(x_2-x_1)/(x_3-x_2)) / 6 f(x_3) + * + * @param[in] idx_range + * The index range on which the quadrature will be carried out. + * + * @return The quadrature coefficients for the Simpson method defined on the provided index range. + */ +template +DFieldMem, typename ExecSpace::memory_space> simpson_quadrature_coefficients_1d( + IdxRange const& idx_range) +{ + DFieldMem, typename ExecSpace::memory_space> coefficients_alloc(idx_range); + fill_simpson_quadrature_coefficients_1d(get_field(coefficients_alloc)); return coefficients_alloc; } +/** + * @brief Get the Simpson coefficients in 1D. + * + * If the number of grid points is not compatible with the Simpson quadrature scheme then + * use a trapezoid formula over one cell at the specified extremity. + * + * @param[in] idx_range + * The index range on which the quadrature will be carried out. + * @param[in] trapezoid_extremity + * The extremity where the trapezoid quadrature may be used. + * + * @return The quadrature coefficients for the Simpson method defined on the provided index range. + */ +template +DFieldMem, typename ExecSpace::memory_space> +simpson_trapezoid_quadrature_coefficients_1d( + IdxRange const& idx_range, + Extremity trapezoid_extremity) +{ + static_assert( + !Grid1D::continuous_dimension_type::PERIODIC, + "The extremity is non-sensical in a Periodic dimension"); + try { + return simpson_quadrature_coefficients_1d(idx_range); + } catch (const std::runtime_error& error) { + DFieldMem, typename ExecSpace::memory_space> coefficients_alloc(idx_range); + DField, typename ExecSpace::memory_space> coefficients( + get_field(coefficients_alloc)); + IdxStep npts_to_remove(1); + if (trapezoid_extremity == Extremity::FRONT) { + fill_simpson_quadrature_coefficients_1d( + coefficients_alloc[idx_range.remove_first(npts_to_remove)]); + Kokkos::parallel_for( + "trapezoid_region", + Kokkos::RangePolicy(0, 1), + KOKKOS_LAMBDA(const int i) { + Idx idx_r = idx_range.back(); + double dx_cell = distance_at_left(idx_r); + coefficients(idx_r) = dx_cell; + coefficients(idx_r - 1) += dx_cell; + }); + } else { + fill_simpson_quadrature_coefficients_1d( + coefficients_alloc[idx_range.remove_last(npts_to_remove)]); + Kokkos::parallel_for( + "trapezoid_region", + Kokkos::RangePolicy(0, 1), + KOKKOS_LAMBDA(const int i) { + Idx idx_l = idx_range.front(); + double dx_cell = distance_at_right(idx_l); + coefficients(idx_l) = dx_cell; + coefficients(idx_l + 1) += dx_cell; + }); + } + return coefficients_alloc; + } +} + /** * @brief Get the simpson coefficients in ND. * diff --git a/src/utils/geometry_descriptors.hpp b/src/utils/geometry_descriptors.hpp new file mode 100644 index 000000000..e1ebde3d5 --- /dev/null +++ b/src/utils/geometry_descriptors.hpp @@ -0,0 +1,4 @@ +// SPDX-License-Identifier: MIT +#pragma once + +enum Extremity { FRONT, BACK };