diff --git a/src/quadrature/CMakeLists.txt b/src/quadrature/CMakeLists.txt index af05e22ac..33d35c674 100644 --- a/src/quadrature/CMakeLists.txt +++ b/src/quadrature/CMakeLists.txt @@ -7,6 +7,7 @@ target_include_directories(quadrature target_link_libraries(quadrature INTERFACE DDC::DDC sll::SLL + gslx::utils ) add_library(gslx::quadrature ALIAS quadrature) diff --git a/src/quadrature/neumann_spline_quadrature.hpp b/src/quadrature/neumann_spline_quadrature.hpp index 57a85b961..56b70cf48 100644 --- a/src/quadrature/neumann_spline_quadrature.hpp +++ b/src/quadrature/neumann_spline_quadrature.hpp @@ -14,11 +14,13 @@ #include +#include "ddc_aliases.hpp" + namespace { -template -using CoefficientChunk1D = ddc::Chunk>; +template +using CoefficientChunk1D = host_t>>; } @@ -33,16 +35,16 @@ using CoefficientChunk1D = ddc::Chunk>; * [1] Non-Uniform Numerical Schemes for the Modelling of Turbulence in the 5D GYSELA Code * Emily Bourne, December 2022- * - * @param[in] domain - * The domain on which the splines quadrature will be carried out. + * @param[in] idx_range + * The index range on which the splines quadrature will be carried out. * @param[in] builder * The spline builder used for the quadrature coefficients. * - * @return The quadrature coefficients for the method defined on the provided domain. + * @return The quadrature coefficients for the method defined on the provided index range. */ -template -ddc::Chunk> neumann_spline_quadrature_coefficients_1d( - ddc::DiscreteDomain const& domain, +template +host_t>> neumann_spline_quadrature_coefficients_1d( + IdxRange const& idx_range, SplineBuilder const& builder) { constexpr int nbc_xmin = SplineBuilder::s_nbc_xmin; @@ -66,12 +68,12 @@ ddc::Chunk> neumann_spline_quadrature_coeffici using bsplines_type = typename SplineBuilder::bsplines_type; - assert(domain.size() == ddc::discrete_space().nbasis() - nbc_xmin - nbc_xmax); + assert(idx_range.size() == ddc::discrete_space().nbasis() - nbc_xmin - nbc_xmax); // Vector of integrals of B-splines - ddc::Chunk> integral_bsplines( - builder.spline_domain()); - ddc::discrete_space().integrals(integral_bsplines.span_view()); + host_t>> integral_bsplines( + get_spline_idx_range(builder)); + ddc::discrete_space().integrals(get_field(integral_bsplines)); // Solve matrix equation Kokkos::View @@ -89,14 +91,13 @@ ddc::Chunk> neumann_spline_quadrature_coeffici .solve(integral_bsplines_mirror_with_additional_allocation, true); Kokkos::deep_copy(integral_bsplines.allocation_kokkos_view(), integral_bsplines_mirror); - ddc::Chunk> coefficients(domain); + host_t>> coefficients(idx_range); // Coefficients of quadrature in integral_bsplines (values which would always be multiplied // by f'(x)=0 are removed - ddc::DiscreteDomain slice - = builder.spline_domain() - .remove(ddc::DiscreteVector {nbc_xmin}, - ddc::DiscreteVector {nbc_xmax}); + IdxRange slice + = get_spline_idx_range(builder) + .remove(IdxStep {nbc_xmin}, IdxStep {nbc_xmax}); Kokkos::deep_copy( coefficients.allocation_kokkos_view(), @@ -114,16 +115,16 @@ ddc::Chunk> neumann_spline_quadrature_coeffici * to calculating and integrating a spline approximation of a function. The spline approximation * would be calculated with homogeneous Neumann boundary conditions. * - * @param[in] domain - * The domain on which the coefficients will be defined. + * @param[in] idx_range + * The index range on which the coefficients will be defined. * @param[in] builders * The spline builder used for the quadrature coefficients in the different dimensions. * * @return The coefficients which define the spline quadrature method in ND. */ template -ddc::Chunk> neumann_spline_quadrature_coefficients( - ddc::DiscreteDomain const& domain, +host_t>> neumann_spline_quadrature_coefficients( + IdxRange const& idx_range, SplineBuilders const&... builders) { assert((std::is_same_v< @@ -132,12 +133,12 @@ ddc::Chunk> neumann_spline_quadrature_coef // Get coefficients for each dimension std::tuple...> current_dim_coeffs( - neumann_spline_quadrature_coefficients_1d(ddc::select(domain), builders)...); + neumann_spline_quadrature_coefficients_1d(ddc::select(idx_range), builders)...); // Allocate ND coefficients - ddc::Chunk> coefficients(domain); + host_t>> coefficients(idx_range); - ddc::for_each(domain, [&](ddc::DiscreteElement const idim) { + ddc::for_each(idx_range, [&](Idx const idim) { // multiply the 1D coefficients by one another coefficients(idim) = (std::get>(current_dim_coeffs)(ddc::select(idim)) diff --git a/src/quadrature/quadrature.hpp b/src/quadrature/quadrature.hpp index c6bbda4e1..d3d659a11 100644 --- a/src/quadrature/quadrature.hpp +++ b/src/quadrature/quadrature.hpp @@ -7,32 +7,33 @@ #include +#include "ddc_aliases.hpp" #include "ddc_helper.hpp" /** - * @brief A class providing an operator for integrating functions defined on a discrete domain. + * @brief A class providing an operator for integrating functions defined on a discrete index range. * - * @tparam QuadratureDomain The domain over which the function is integrated. - * @tparam TotalDomain The domain of the chunk which can be passed to the operator(). This is the - * QuadratureDomain combined with any batch dimensions. If there are no + * @tparam QuadratureIdxRange The index range over which the function is integrated. + * @tparam TotalIdxRange The index range of the chunk which can be passed to the operator(). This is the + * QuadratureIdxRange combined with any batch dimensions. If there are no * batch dimensions then this argument does not need to be provided as by - * default it is equal to the QuadratureDomain. + * default it is equal to the QuadratureIdxRange. * @tparam MemorySpace The memory space (cpu/gpu) where the quadrature coefficients are saved. */ template < - class QuadratureDomain, - class TotalDomain = QuadratureDomain, + class QuadratureIdxRange, + class TotalIdxRange = QuadratureIdxRange, class MemorySpace = Kokkos::DefaultExecutionSpace::memory_space> class Quadrature { private: /// The tyoe of an element of an index of the quadrature coefficients. - using QuadratureIndex = typename QuadratureDomain::discrete_element_type; + using QuadratureIdx = typename QuadratureIdxRange::discrete_element_type; - using QuadChunkView = ddc:: - ChunkView; + using QuadConstField + = DConstField; - QuadChunkView m_coefficients; + QuadConstField m_coefficients; public: /** @@ -40,20 +41,20 @@ class Quadrature * @param coeffs * The coefficients of the quadrature. */ - explicit Quadrature(QuadChunkView coeffs) : m_coefficients(coeffs) {} + explicit Quadrature(QuadConstField coeffs) : m_coefficients(coeffs) {} /** - * @brief An operator for calculating the integral of a function defined on a discrete domain. + * @brief An operator for calculating the integral of a function defined on a discrete index range. * * @param[in] exec_space * The space on which the function is executed (CPU/GPU). * @param[in] integrated_function - * A function taking an index of a position in the domain over which the quadrature is + * A function taking an index of a position in the index range over which the quadrature is * calculated and returning the value of the function to be integrated at that point. * It should be noted that a ChunkSpan fulfils these criteria and can be passed as the function to be integrated. * If the exec_space is a GPU the function that is passed must be accessible from GPU. * - * @returns The integral of the function over the domain. + * @returns The integral of the function over the index range. */ template double operator()(ExecutionSpace exec_space, IntegratorFunction integrated_function) const @@ -62,9 +63,9 @@ class Quadrature Kokkos::SpaceAccessibility::accessible, "Execution space is not compatible with memory space where coefficients are found"); static_assert( - std::is_invocable_v, + std::is_invocable_v, "The object passed to Quadrature::operator() is not defined on the quadrature " - "domain."); + "idx_range."); ddc::ChunkSpan const coeff_proxy = m_coefficients; @@ -72,16 +73,16 @@ class Quadrature exec_space.fence(); return ddc::parallel_transform_reduce( exec_space, - coeff_proxy.domain(), + get_idx_range(coeff_proxy), 0.0, ddc::reducer::sum(), - KOKKOS_LAMBDA(QuadratureIndex const ix) { + KOKKOS_LAMBDA(QuadratureIdx const ix) { return coeff_proxy(ix) * integrated_function(ix); }); } /** - * @brief An operator for calculating the integral of a function defined on a discrete domain + * @brief An operator for calculating the integral of a function defined on a discrete index range * by cycling over batch dimensions. * * @param[in] exec_space @@ -89,17 +90,16 @@ class Quadrature * @param[out] result * The result of the quadrature calculation. * @param[in] integrated_function - * A function taking an index of a position in the domain over which the quadrature is - * calculated (including the batch domain) and returning the value of the function to + * A function taking an index of a position in the index range over which the quadrature is + * calculated (including the batch index range) and returning the value of the function to * be integrated at that point. * Please note that a ChunkSpan fulfills the described criteria. * If the exec_space is a GPU the function that is passed must be accessible from GPU. */ - template + template void operator()( ExecutionSpace exec_space, - ddc::ChunkSpan const - result, + Field const result, IntegratorFunction integrated_function) const { static_assert( @@ -110,40 +110,42 @@ class Quadrature "Kokkos::TeamPolicy only works with the default execution space. Please use " "DefaultExecutionSpace to call this batched operator."); using ExpectedBatchDims = ddc::type_seq_remove_t< - ddc::to_type_seq_t, - ddc::to_type_seq_t>; + ddc::to_type_seq_t, + ddc::to_type_seq_t>; static_assert( - ddc::type_seq_same_v, ExpectedBatchDims>, - "The batch domain deduced from the type of result does not match the class " + ddc::type_seq_same_v, ExpectedBatchDims>, + "The batch idx_range deduced from the type of result does not match the class " "template parameters."); // Get useful index types - using TotalIndex = typename TotalDomain::discrete_element_type; - using BatchIndex = typename BatchDomain::discrete_element_type; + using TotalIdx = typename TotalIdxRange::discrete_element_type; + using BatchIdx = typename BatchIdxRange::discrete_element_type; static_assert( - std::is_invocable_v, - "The object passed to Quadrature::operator() is not defined on the total domain."); + std::is_invocable_v, + "The object passed to Quadrature::operator() is not defined on the total " + "idx_range."); - // Get domains - QuadratureDomain quad_domain(m_coefficients.domain()); - BatchDomain batch_domain(result.domain()); + // Get index ranges + QuadratureIdxRange quad_idx_range(get_idx_range(m_coefficients)); + BatchIdxRange batch_idx_range(get_idx_range(result)); ddc::ChunkSpan const coeff_proxy = m_coefficients; // Loop over batch dimensions Kokkos::parallel_for( - Kokkos::TeamPolicy<>(exec_space, batch_domain.size(), Kokkos::AUTO), + Kokkos::TeamPolicy<>(exec_space, batch_idx_range.size(), Kokkos::AUTO), KOKKOS_LAMBDA(const Kokkos::TeamPolicy<>::member_type& team) { const int idx = team.league_rank(); - BatchIndex ib = to_discrete_element(idx, batch_domain); + BatchIdx ib = to_discrete_element(idx, batch_idx_range); // Sum over quadrature dimensions double teamSum = 0; Kokkos::parallel_reduce( - Kokkos::TeamThreadRange(team, quad_domain.size()), + Kokkos::TeamThreadRange(team, quad_idx_range.size()), [&](int const& thread_index, double& sum) { - QuadratureIndex iq = to_discrete_element(thread_index, quad_domain); - TotalIndex it(ib, iq); + QuadratureIdx iq + = to_discrete_element(thread_index, quad_idx_range); + TotalIdx it(ib, iq); sum += coeff_proxy(iq) * integrated_function(it); }, teamSum); @@ -153,37 +155,35 @@ class Quadrature private: /** - * A function which converts an integer into a DiscreteElement found in a domain - * starting from the front. This is useful for iterating over a domain using Kokkos + * A function which converts an integer into a DiscreteElement found in an index range + * starting from the front. This is useful for iterating over an index range using Kokkos * loops. * * @param[in] idx The index of the requested element. - * @param[in] dom The domain being iterated over. + * @param[in] dom The index range being iterated over. * - * @return The vector displacement from the front of the domain + * @return The vector displacement from the front of the index range */ template - KOKKOS_FUNCTION static ddc::DiscreteElement to_discrete_element( + KOKKOS_FUNCTION static Idx to_discrete_element( int idx, - ddc::DiscreteDomain dom) + IdxRange dom) { - ddc::DiscreteDomain subdomain(dom); - ddc::DiscreteElement head_idx( - ddc::select(dom).front() + idx / subdomain.size()); + IdxRange subidx_range(dom); + Idx head_idx(ddc::select(dom).front() + idx / subidx_range.size()); if constexpr (sizeof...(Grid1D) == 0) { return head_idx; } else { - ddc::DiscreteElement tail_idx - = to_discrete_element(idx % subdomain.size(), subdomain); - return ddc::DiscreteElement(head_idx, tail_idx); + Idx tail_idx = to_discrete_element(idx % subidx_range.size(), subidx_range); + return Idx(head_idx, tail_idx); } } }; namespace detail { -template -struct OnMemorySpace> +template +struct OnMemorySpace> { - using type = Quadrature; + using type = Quadrature; }; } // namespace detail diff --git a/src/quadrature/quadrature_coeffs_nd.hpp b/src/quadrature/quadrature_coeffs_nd.hpp index e76993b21..4b8b9f4aa 100644 --- a/src/quadrature/quadrature_coeffs_nd.hpp +++ b/src/quadrature/quadrature_coeffs_nd.hpp @@ -7,36 +7,37 @@ #include +#include "ddc_aliases.hpp" + namespace { -template -using CoefficientChunk1D = ddc::Chunk>; +template +using CoefficientChunk1D = host_t>>; } /** * @brief Helper function which creates ND dimensions from N 1D quadrature coefficient functions. * - * @param domain - * The domain on which the coefficients will be defined. + * @param idx_range + * The index range on which the coefficients will be defined. * @param funcs * The functions which define quadrature coefficients in the different dimensions. * * @returns The coefficients which define the quadrature method in ND. */ template -ddc::Chunk> quadrature_coeffs_nd( - ddc::DiscreteDomain const& domain, - std::function>( - ddc::DiscreteDomain)>... funcs) +host_t>> quadrature_coeffs_nd( + IdxRange const& idx_range, + std::function>>(IdxRange)>... funcs) { // Get coefficients for each dimension std::tuple...> current_dim_coeffs( - funcs(ddc::select(domain))...); + funcs(ddc::select(idx_range))...); // Allocate ND coefficients - ddc::Chunk> coefficients(domain); + host_t>> coefficients(idx_range); - ddc::for_each(domain, [&](ddc::DiscreteElement const idim) { + ddc::for_each(idx_range, [&](Idx const idim) { // multiply the 1D coefficients by one another coefficients(idim) = (std::get>(current_dim_coeffs)(ddc::select(idim)) diff --git a/src/quadrature/simpson_quadrature.hpp b/src/quadrature/simpson_quadrature.hpp index dc65777cc..a127521fb 100644 --- a/src/quadrature/simpson_quadrature.hpp +++ b/src/quadrature/simpson_quadrature.hpp @@ -5,35 +5,37 @@ #include +#include "ddc_aliases.hpp" + /** * @brief Get the Simpson coefficients in 1D. * - * Calculate the quadrature coefficients for the Simpson method defined on the provided domain. + * Calculate the quadrature coefficients for the Simpson method defined on the provided index range. * - * @param[in] domain - * The domain on which the quadrature will be carried out. + * @param[in] index range + * The index range on which the quadrature will be carried out. * - * @return The quadrature coefficients for the Simpson method defined on the provided domain. + * @return The quadrature coefficients for the Simpson method defined on the provided index range. */ -template -ddc::Chunk> simpson_quadrature_coefficients_1d( - ddc::DiscreteDomain const& domain) +template +host_t>> simpson_quadrature_coefficients_1d( + IdxRange const& idx_range) { - ddc::Chunk> coefficients(domain); + host_t>> coefficients(idx_range); - coefficients(domain.front()) = 1. / 3. * distance_at_right(domain.front()); + coefficients(idx_range.front()) = 1. / 3. * distance_at_right(idx_range.front()); - for (auto it = domain.begin() + 1; it < domain.end() - 1; it += 2) { - ddc::DiscreteElement idx = *it; + for (auto it = idx_range.begin() + 1; it < idx_range.end() - 1; it += 2) { + Idx idx = *it; coefficients(idx) = 2. / 3. * (distance_at_left(idx) + distance_at_right(idx)); - idx += ddc::DiscreteVector(1); + idx += IdxStep(1); coefficients(idx) = 1. / 3. * (distance_at_left(idx) + distance_at_right(idx)); } - coefficients(domain.back()) = 1. / 3. * distance_at_left(domain.back()); + coefficients(idx_range.back()) = 1. / 3. * distance_at_left(idx_range.back()); - if constexpr (IDim::continuous_dimension_type::PERIODIC) { - coefficients(domain.front()) += 2. / 3 * distance_at_left(domain.back()); - coefficients(domain.back()) += 2. / 3 * distance_at_right(domain.front()); + if constexpr (Grid::continuous_dimension_type::PERIODIC) { + coefficients(idx_range.front()) += 2. / 3 * distance_at_left(idx_range.back()); + coefficients(idx_range.back()) += 2. / 3 * distance_at_right(idx_range.front()); } return coefficients; diff --git a/src/quadrature/spline_quadrature.hpp b/src/quadrature/spline_quadrature.hpp index e35dcb99e..9d46f8954 100644 --- a/src/quadrature/spline_quadrature.hpp +++ b/src/quadrature/spline_quadrature.hpp @@ -13,11 +13,13 @@ #include +#include "ddc_aliases.hpp" + namespace { -template -using CoefficientChunk1D = ddc::Chunk>; +template +using CoefficientChunk1D = host_t>>; } @@ -46,17 +48,17 @@ using CoefficientChunk1D = ddc::Chunk>; * in the 5D GYSELA Code". December 2022. * * - * @param[in] domain - * The domain where the functions we want to integrate + * @param[in] idx_range + * The index range where the functions we want to integrate * are defined. * @param[in] builder * The spline builder describing the way in which splines would be constructed. * * @return A chunk with the quadrature coefficients @f$ q @f$. */ -template -ddc::Chunk> spline_quadrature_coefficients_1d( - ddc::DiscreteDomain const& domain, +template +host_t>> spline_quadrature_coefficients_1d( + IdxRange const& idx_range, SplineBuilder const& builder) { static_assert( @@ -71,15 +73,15 @@ ddc::Chunk> spline_quadrature_coefficients_1d( using bsplines_type = typename SplineBuilder::bsplines_type; // Vector of integrals of B-splines - ddc::Chunk> integral_bsplines( - builder.spline_domain()); - ddc::discrete_space().integrals(integral_bsplines.span_view()); + host_t>> integral_bsplines( + get_spline_idx_range(builder)); + ddc::discrete_space().integrals(get_field(integral_bsplines)); // Solve matrix equation ddc::ChunkSpan integral_bsplines_without_periodic_point - = integral_bsplines.span_view()[ddc::DiscreteDomain( - ddc::DiscreteElement(0), - ddc::DiscreteVector(builder.get_interpolation_matrix().size()))]; + = integral_bsplines[IdxRange( + Idx(0), + IdxStep(builder.get_interpolation_matrix().size()))]; Kokkos::View integral_bsplines_mirror_with_additional_allocation( "integral_bsplines_mirror_with_additional_allocation", @@ -101,7 +103,7 @@ ddc::Chunk> spline_quadrature_coefficients_1d( integral_bsplines_without_periodic_point.allocation_kokkos_view(), integral_bsplines_mirror); - ddc::Chunk> coefficients(domain); + host_t>> coefficients(idx_range); Kokkos::deep_copy( coefficients.allocation_kokkos_view(), @@ -115,18 +117,18 @@ ddc::Chunk> spline_quadrature_coefficients_1d( /** * @brief Get the spline quadrature coefficients in ND from N 1D quadrature coefficient. * - * Calculate the quadrature coefficients for the spline quadrature method defined on the provided domain. + * Calculate the quadrature coefficients for the spline quadrature method defined on the provided index range. * - * @param[in] domain - * The domain on which the coefficients will be defined. + * @param[in] idx_range + * The index range on which the coefficients will be defined. * @param[in] builders * The spline builder used for the quadrature coefficients in the different dimensions. * * @return The coefficients which define the spline quadrature method in ND. */ template -ddc::Chunk> spline_quadrature_coefficients( - ddc::DiscreteDomain const& domain, +host_t>> spline_quadrature_coefficients( + IdxRange const& idx_range, SplineBuilders const&... builders) { assert((std::is_same_v< @@ -135,12 +137,12 @@ ddc::Chunk> spline_quadrature_coefficients // Get coefficients for each dimension std::tuple...> current_dim_coeffs( - spline_quadrature_coefficients_1d(ddc::select(domain), builders)...); + spline_quadrature_coefficients_1d(ddc::select(idx_range), builders)...); // Allocate ND coefficients - ddc::Chunk> coefficients(domain); + host_t>> coefficients(idx_range); - ddc::for_each(domain, [&](ddc::DiscreteElement const idim) { + ddc::for_each(idx_range, [&](Idx const idim) { // multiply the 1D coefficients by one another coefficients(idim) = (std::get>(current_dim_coeffs)(ddc::select(idim)) diff --git a/src/quadrature/trapezoid_quadrature.hpp b/src/quadrature/trapezoid_quadrature.hpp index c38d27797..397fb60a3 100644 --- a/src/quadrature/trapezoid_quadrature.hpp +++ b/src/quadrature/trapezoid_quadrature.hpp @@ -3,38 +3,37 @@ * File providing quadrature coefficients via the trapezoidal method. */ #pragma once - #include +#include "ddc_aliases.hpp" #include "quadrature_coeffs_nd.hpp" /** * @brief Get the trapezoid coefficients in 1D. * - * Calculate the quadrature coefficients for the trapezoid method defined on the provided domain. + * Calculate the quadrature coefficients for the trapezoid method defined on the provided index range. * - * @param[in] domain - * The domain on which the quadrature will be carried out. + * @param[in] idx_range + * The index range on which the quadrature will be carried out. * - * @return The quadrature coefficients for the trapezoid method defined on the provided domain. + * @return The quadrature coefficients for the trapezoid method defined on the provided index range. */ -template -ddc::Chunk> trapezoid_quadrature_coefficients_1d( - ddc::DiscreteDomain const& domain) +template +host_t>> trapezoid_quadrature_coefficients_1d( + IdxRange const& idx_range) { - ddc::Chunk> coefficients(domain); - ddc::DiscreteDomain middle_domain - = domain.remove(ddc::DiscreteVector(1), ddc::DiscreteVector(1)); + host_t>> coefficients(idx_range); + IdxRange middle_idx_range = idx_range.remove(IdxStep(1), IdxStep(1)); - coefficients(domain.front()) = 0.5 * distance_at_right(domain.front()); - ddc::for_each(middle_domain, [&](ddc::DiscreteElement const idx) { + coefficients(idx_range.front()) = 0.5 * distance_at_right(idx_range.front()); + ddc::for_each(middle_idx_range, [&](Idx const idx) { coefficients(idx) = 0.5 * (distance_at_left(idx) + distance_at_right(idx)); }); - coefficients(domain.back()) = 0.5 * distance_at_left(domain.back()); + coefficients(idx_range.back()) = 0.5 * distance_at_left(idx_range.back()); - if constexpr (IDim::continuous_dimension_type::PERIODIC) { - coefficients(domain.front()) += 0.5 * distance_at_left(domain.back()); - coefficients(domain.back()) += 0.5 * distance_at_right(domain.front()); + if constexpr (Grid::continuous_dimension_type::PERIODIC) { + coefficients(idx_range.front()) += 0.5 * distance_at_left(idx_range.back()); + coefficients(idx_range.back()) += 0.5 * distance_at_right(idx_range.front()); } return coefficients; @@ -43,19 +42,19 @@ ddc::Chunk> trapezoid_quadrature_coefficients_ /** * @brief Get the trapezoid coefficients in ND. * - * Calculate the quadrature coefficients for the trapezoid method defined on the provided domain. + * Calculate the quadrature coefficients for the trapezoid method defined on the provided index range. * - * @param[in] domain - * The domain on which the quadrature will be carried out. + * @param[in] idx_range + * The index range on which the quadrature will be carried out. * - * @return The quadrature coefficients for the trapezoid method defined on the provided domain. + * @return The quadrature coefficients for the trapezoid method defined on the provided index range. */ template -ddc::Chunk> trapezoid_quadrature_coefficients( - ddc::DiscreteDomain const& domain) +host_t>> trapezoid_quadrature_coefficients( + IdxRange const& idx_range) { return quadrature_coeffs_nd( - domain, - (std::function>( - ddc::DiscreteDomain)>(trapezoid_quadrature_coefficients_1d))...); + idx_range, + (std::function>>(IdxRange)>( + trapezoid_quadrature_coefficients_1d))...); } diff --git a/tests/quadrature/batched_quadrature.cpp b/tests/quadrature/batched_quadrature.cpp index 0cb4d7edc..4c5a8dec8 100644 --- a/tests/quadrature/batched_quadrature.cpp +++ b/tests/quadrature/batched_quadrature.cpp @@ -29,94 +29,94 @@ struct Y static bool constexpr PERIODIC = false; }; -using CoordBatch1 = const ddc::Coordinate; -using CoordBatch2 = const ddc::Coordinate; -using CoordX = const ddc::Coordinate; -using CoordY = const ddc::Coordinate; +using CoordBatch1 = const Coord; +using CoordBatch2 = const Coord; +using CoordX = const Coord; +using CoordY = const Coord; -struct IDimBatch1 : ddc::UniformPointSampling +struct GridBatch1 : UniformGridBase { }; -struct IDimBatch2 : ddc::UniformPointSampling +struct GridBatch2 : UniformGridBase { }; -struct IDimX : ddc::UniformPointSampling +struct GridX : UniformGridBase { }; -struct IDimY : ddc::UniformPointSampling +struct GridY : UniformGridBase { }; -using IdxBatch1 = ddc::DiscreteElement; -using IdxBatch2 = ddc::DiscreteElement; -using IdxX = ddc::DiscreteElement; -using IdxY = ddc::DiscreteElement; -using IdxXY = ddc::DiscreteElement; -using IdxB1X = ddc::DiscreteElement; -using IdxB1B2 = ddc::DiscreteElement; -using IdxB1XY = ddc::DiscreteElement; -using IdxB1B2X = ddc::DiscreteElement; -using IdxB1B2XY = ddc::DiscreteElement; -using IdxXB1YB2 = ddc::DiscreteElement; - -using IVectBatch1 = const ddc::DiscreteVector; -using IVectBatch2 = const ddc::DiscreteVector; -using IVectX = const ddc::DiscreteVector; -using IVectY = const ddc::DiscreteVector; - -using IDomainBatch1 = ddc::DiscreteDomain; -using IDomainBatch2 = ddc::DiscreteDomain; -using IDomainX = ddc::DiscreteDomain; -using IDomainY = ddc::DiscreteDomain; -using IDomainB1X = ddc::DiscreteDomain; -using IDomainB1XY = ddc::DiscreteDomain; -using IDomainB1B2X = ddc::DiscreteDomain; -using IDomainB1B2XY = ddc::DiscreteDomain; -using IDomainXB1YB2 = ddc::DiscreteDomain; -using IDomainXY = ddc::DiscreteDomain; -using IDomainB1B2 = ddc::DiscreteDomain; - -using DFieldBatch1 = device_t>; -using DFieldX = device_t>; -using DFieldY = device_t>; -using DFieldXY = device_t>; -using DFieldB1B2 = device_t>; +using IdxBatch1 = Idx; +using IdxBatch2 = Idx; +using IdxX = Idx; +using IdxY = Idx; +using IdxXY = Idx; +using IdxB1X = Idx; +using IdxB1B2 = Idx; +using IdxB1XY = Idx; +using IdxB1B2X = Idx; +using IdxB1B2XY = Idx; +using IdxXB1YB2 = Idx; + +using IdxStepBatch1 = const IdxStep; +using IdxStepBatch2 = const IdxStep; +using IdxStepX = const IdxStep; +using IdxStepY = const IdxStep; + +using IdxRangeBatch1 = IdxRange; +using IdxRangeBatch2 = IdxRange; +using IdxRangeX = IdxRange; +using IdxRangeY = IdxRange; +using IdxRangeB1X = IdxRange; +using IdxRangeB1XY = IdxRange; +using IdxRangeB1B2X = IdxRange; +using IdxRangeB1B2XY = IdxRange; +using IdxRangeXB1YB2 = IdxRange; +using IdxRangeXY = IdxRange; +using IdxRangeB1B2 = IdxRange; + +using DFieldMemBatch1 = FieldMem; +using DFieldMemX = FieldMem; +using DFieldMemY = FieldMem; +using DFieldMemXY = FieldMem; +using DFieldMemB1B2 = FieldMem; void batched_operator_1d() { CoordBatch1 b_min(0.0); CoordBatch1 b_max(3.0); - IVectBatch1 b_ncells(4); + IdxStepBatch1 b_ncells(4); CoordX x_min(4.0); CoordX x_max(8.0); - IVectX x_ncells(16); + IdxStepX x_ncells(16); - IDomainBatch1 gridb = ddc::init_discrete_space( - IDimBatch1::init(b_min, b_max, b_ncells)); - IDomainX gridx = ddc::init_discrete_space(IDimX::init(x_min, x_max, x_ncells)); + IdxRangeBatch1 gridb = ddc::init_discrete_space( + GridBatch1::init(b_min, b_max, b_ncells)); + IdxRangeX gridx = ddc::init_discrete_space(GridX::init(x_min, x_max, x_ncells)); - host_t quad_coeffs_host = trapezoid_quadrature_coefficients(gridx); + host_t quad_coeffs_host = trapezoid_quadrature_coefficients(gridx); auto quad_coeffs = ddc::create_mirror_view_and_copy( Kokkos::DefaultExecutionSpace(), - quad_coeffs_host.span_view()); - Quadrature quad_batched_operator(quad_coeffs.span_view()); + get_field(quad_coeffs_host)); + Quadrature quad_batched_operator(get_field(quad_coeffs)); - DFieldBatch1 results(gridb); + DFieldMemBatch1 results(gridb); quad_batched_operator( Kokkos::DefaultExecutionSpace(), - results.span_view(), + get_field(results), KOKKOS_LAMBDA(IdxB1X ibx) { - double b = ddc::coordinate(ddc::select(ibx)); - double x = ddc::coordinate(ddc::select(ibx)); + double b = ddc::coordinate(ddc::select(ibx)); + double x = ddc::coordinate(ddc::select(ibx)); return b * x + 2; }); auto results_host = ddc:: - create_mirror_view_and_copy(Kokkos::DefaultHostExecutionSpace(), results.span_view()); + create_mirror_view_and_copy(Kokkos::DefaultHostExecutionSpace(), get_field(results)); ddc::for_each(gridb, [&](IdxBatch1 ib) { - double b = ddc::coordinate(ddc::select(ib)); + double b = ddc::coordinate(ddc::select(ib)); double x = x_max; double const ubound = 0.5 * b * x * x + 2 * x; x = x_min; @@ -129,52 +129,52 @@ void batched_operator_2d() { CoordBatch1 b1_min(0.0); CoordBatch1 b1_max(3.0); - IVectBatch1 b1_ncells(4); + IdxStepBatch1 b1_ncells(4); CoordBatch2 b2_min(1.0); CoordBatch2 b2_max(2.0); - IVectBatch2 b2_ncells(3); + IdxStepBatch2 b2_ncells(3); CoordX x_min(4.0); CoordX x_max(8.0); - IVectX x_ncells(16); + IdxStepX x_ncells(16); CoordY y_min(4.0); CoordY y_max(8.0); - IVectY y_ncells(16); + IdxStepY y_ncells(16); - IDomainBatch1 gridb1 = ddc::init_discrete_space( - IDimBatch1::init(b1_min, b1_max, b1_ncells)); - IDomainBatch2 gridb2 = ddc::init_discrete_space( - IDimBatch2::init(b2_min, b2_max, b2_ncells)); - IDomainX gridx = ddc::init_discrete_space(IDimX::init(x_min, x_max, x_ncells)); - IDomainY gridy = ddc::init_discrete_space(IDimY::init(y_min, y_max, y_ncells)); + IdxRangeBatch1 gridb1 = ddc::init_discrete_space( + GridBatch1::init(b1_min, b1_max, b1_ncells)); + IdxRangeBatch2 gridb2 = ddc::init_discrete_space( + GridBatch2::init(b2_min, b2_max, b2_ncells)); + IdxRangeX gridx = ddc::init_discrete_space(GridX::init(x_min, x_max, x_ncells)); + IdxRangeY gridy = ddc::init_discrete_space(GridY::init(y_min, y_max, y_ncells)); - IDomainXY gridxy(gridx, gridy); + IdxRangeXY gridxy(gridx, gridy); - host_t quad_coeffs_host = trapezoid_quadrature_coefficients(gridxy); + host_t quad_coeffs_host = trapezoid_quadrature_coefficients(gridxy); auto quad_coeffs = ddc::create_mirror_view_and_copy( Kokkos::DefaultExecutionSpace(), - quad_coeffs_host.span_view()); - Quadrature quad_batched_operator(quad_coeffs.span_view()); + get_field(quad_coeffs_host)); + Quadrature quad_batched_operator(get_field(quad_coeffs)); - IDomainB1B2 gridb(gridb1, gridb2); + IdxRangeB1B2 gridb(gridb1, gridb2); - DFieldB1B2 results(gridb); + DFieldMemB1B2 results(gridb); quad_batched_operator( Kokkos::DefaultExecutionSpace(), - results.span_view(), + get_field(results), KOKKOS_LAMBDA(IdxB1B2XY ibx) { - double const x = ddc::coordinate(ddc::select(ibx)); - double const y = ddc::coordinate(ddc::select(ibx)); - double const b1 = ddc::coordinate(ddc::select(ibx)); - double const b2 = ddc::coordinate(ddc::select(ibx)); + double const x = ddc::coordinate(ddc::select(ibx)); + double const y = ddc::coordinate(ddc::select(ibx)); + double const b1 = ddc::coordinate(ddc::select(ibx)); + double const b2 = ddc::coordinate(ddc::select(ibx)); return b1 * (x * y + 2 * x + b2 * y); }); auto results_host = ddc:: - create_mirror_view_and_copy(Kokkos::DefaultHostExecutionSpace(), results.span_view()); + create_mirror_view_and_copy(Kokkos::DefaultHostExecutionSpace(), get_field(results)); ddc::for_each(gridb, [&](IdxB1B2 ib) { - double const b1 = ddc::coordinate(ddc::select(ib)); - double const b2 = ddc::coordinate(ddc::select(ib)); + double const b1 = ddc::coordinate(ddc::select(ib)); + double const b2 = ddc::coordinate(ddc::select(ib)); double const exact = b1 * (y_max * y_max @@ -192,40 +192,40 @@ void batched_operator_1d_2d() { CoordBatch1 b1_min(0.0); CoordBatch1 b1_max(3.0); - IVectBatch1 b1_ncells(4); + IdxStepBatch1 b1_ncells(4); CoordX x_min(4.0); CoordX x_max(8.0); - IVectX x_ncells(16); + IdxStepX x_ncells(16); CoordY y_min(4.0); CoordY y_max(8.0); - IVectY y_ncells(16); + IdxStepY y_ncells(16); - IDomainBatch1 gridb1 = ddc::init_discrete_space( - IDimBatch1::init(b1_min, b1_max, b1_ncells)); - IDomainX gridx = ddc::init_discrete_space(IDimX::init(x_min, x_max, x_ncells)); - IDomainY gridy = ddc::init_discrete_space(IDimY::init(y_min, y_max, y_ncells)); + IdxRangeBatch1 gridb1 = ddc::init_discrete_space( + GridBatch1::init(b1_min, b1_max, b1_ncells)); + IdxRangeX gridx = ddc::init_discrete_space(GridX::init(x_min, x_max, x_ncells)); + IdxRangeY gridy = ddc::init_discrete_space(GridY::init(y_min, y_max, y_ncells)); - IDomainXY gridxy(gridx, gridy); + IdxRangeXY gridxy(gridx, gridy); - host_t quad_coeffs_host = trapezoid_quadrature_coefficients(gridxy); + host_t quad_coeffs_host = trapezoid_quadrature_coefficients(gridxy); auto quad_coeffs = ddc::create_mirror_view_and_copy( Kokkos::DefaultExecutionSpace(), - quad_coeffs_host.span_view()); - Quadrature quad_batched_operator(quad_coeffs.span_view()); + get_field(quad_coeffs_host)); + Quadrature quad_batched_operator(get_field(quad_coeffs)); - DFieldBatch1 results(gridb1); + DFieldMemBatch1 results(gridb1); quad_batched_operator( Kokkos::DefaultExecutionSpace(), - results.span_view(), + get_field(results), KOKKOS_LAMBDA(IdxB1XY ibx) { - double const x = ddc::coordinate(ddc::select(ibx)); - double const y = ddc::coordinate(ddc::select(ibx)); - double const b1 = ddc::coordinate(ddc::select(ibx)); + double const x = ddc::coordinate(ddc::select(ibx)); + double const y = ddc::coordinate(ddc::select(ibx)); + double const b1 = ddc::coordinate(ddc::select(ibx)); return b1 * (x * y + 2 * x + 3 * y); }); auto results_host = ddc:: - create_mirror_view_and_copy(Kokkos::DefaultHostExecutionSpace(), results.span_view()); + create_mirror_view_and_copy(Kokkos::DefaultHostExecutionSpace(), get_field(results)); ddc::for_each(gridb1, [&](IdxBatch1 ib) { double const b1 = ddc::coordinate(ib); @@ -244,45 +244,45 @@ void batched_operator_2d_1d() { CoordBatch1 b1_min(0.0); CoordBatch1 b1_max(3.0); - IVectBatch1 b1_ncells(4); + IdxStepBatch1 b1_ncells(4); CoordBatch2 b2_min(1.0); CoordBatch2 b2_max(2.0); - IVectBatch2 b2_ncells(3); + IdxStepBatch2 b2_ncells(3); CoordX x_min(4.0); CoordX x_max(8.0); - IVectX x_ncells(16); + IdxStepX x_ncells(16); - IDomainBatch1 gridb1 = ddc::init_discrete_space( - IDimBatch1::init(b1_min, b1_max, b1_ncells)); - IDomainBatch2 gridb2 = ddc::init_discrete_space( - IDimBatch2::init(b2_min, b2_max, b2_ncells)); - IDomainX gridx = ddc::init_discrete_space(IDimX::init(x_min, x_max, x_ncells)); + IdxRangeBatch1 gridb1 = ddc::init_discrete_space( + GridBatch1::init(b1_min, b1_max, b1_ncells)); + IdxRangeBatch2 gridb2 = ddc::init_discrete_space( + GridBatch2::init(b2_min, b2_max, b2_ncells)); + IdxRangeX gridx = ddc::init_discrete_space(GridX::init(x_min, x_max, x_ncells)); - host_t quad_coeffs_host = trapezoid_quadrature_coefficients(gridx); + host_t quad_coeffs_host = trapezoid_quadrature_coefficients(gridx); auto quad_coeffs = ddc::create_mirror_view_and_copy( Kokkos::DefaultExecutionSpace(), - quad_coeffs_host.span_view()); - Quadrature quad_batched_operator(quad_coeffs.span_view()); + get_field(quad_coeffs_host)); + Quadrature quad_batched_operator(get_field(quad_coeffs)); - IDomainB1B2 gridb(gridb1, gridb2); + IdxRangeB1B2 gridb(gridb1, gridb2); - DFieldB1B2 results(gridb); + DFieldMemB1B2 results(gridb); quad_batched_operator( Kokkos::DefaultExecutionSpace(), - results.span_view(), + get_field(results), KOKKOS_LAMBDA(IdxB1B2X ibx) { - double const x = ddc::coordinate(ddc::select(ibx)); - double const b1 = ddc::coordinate(ddc::select(ibx)); - double const b2 = ddc::coordinate(ddc::select(ibx)); + double const x = ddc::coordinate(ddc::select(ibx)); + double const b1 = ddc::coordinate(ddc::select(ibx)); + double const b2 = ddc::coordinate(ddc::select(ibx)); return b1 * x + b2; }); auto results_host = ddc:: - create_mirror_view_and_copy(Kokkos::DefaultHostExecutionSpace(), results.span_view()); + create_mirror_view_and_copy(Kokkos::DefaultHostExecutionSpace(), get_field(results)); ddc::for_each(gridb, [&](IdxB1B2 ib) { - double const b1 = ddc::coordinate(ddc::select(ib)); - double const b2 = ddc::coordinate(ddc::select(ib)); + double const b1 = ddc::coordinate(ddc::select(ib)); + double const b2 = ddc::coordinate(ddc::select(ib)); double x = x_max; double const ubound = 0.5 * b1 * x * x + b2 * x; x = x_min; @@ -295,52 +295,52 @@ void batched_operator_2d_reordered() { CoordBatch1 b1_min(0.0); CoordBatch1 b1_max(3.0); - IVectBatch1 b1_ncells(4); + IdxStepBatch1 b1_ncells(4); CoordBatch2 b2_min(1.0); CoordBatch2 b2_max(2.0); - IVectBatch2 b2_ncells(3); + IdxStepBatch2 b2_ncells(3); CoordX x_min(4.0); CoordX x_max(8.0); - IVectX x_ncells(16); + IdxStepX x_ncells(16); CoordY y_min(4.0); CoordY y_max(8.0); - IVectY y_ncells(16); + IdxStepY y_ncells(16); - IDomainBatch1 gridb1 = ddc::init_discrete_space( - IDimBatch1::init(b1_min, b1_max, b1_ncells)); - IDomainBatch2 gridb2 = ddc::init_discrete_space( - IDimBatch2::init(b2_min, b2_max, b2_ncells)); - IDomainX gridx = ddc::init_discrete_space(IDimX::init(x_min, x_max, x_ncells)); - IDomainY gridy = ddc::init_discrete_space(IDimY::init(y_min, y_max, y_ncells)); + IdxRangeBatch1 gridb1 = ddc::init_discrete_space( + GridBatch1::init(b1_min, b1_max, b1_ncells)); + IdxRangeBatch2 gridb2 = ddc::init_discrete_space( + GridBatch2::init(b2_min, b2_max, b2_ncells)); + IdxRangeX gridx = ddc::init_discrete_space(GridX::init(x_min, x_max, x_ncells)); + IdxRangeY gridy = ddc::init_discrete_space(GridY::init(y_min, y_max, y_ncells)); - IDomainXY gridxy(gridx, gridy); + IdxRangeXY gridxy(gridx, gridy); - host_t quad_coeffs_host = trapezoid_quadrature_coefficients(gridxy); + host_t quad_coeffs_host = trapezoid_quadrature_coefficients(gridxy); auto quad_coeffs = ddc::create_mirror_view_and_copy( Kokkos::DefaultExecutionSpace(), - quad_coeffs_host.span_view()); - Quadrature quad_batched_operator(quad_coeffs.span_view()); + get_field(quad_coeffs_host)); + Quadrature quad_batched_operator(get_field(quad_coeffs)); - IDomainB1B2 gridb(gridb1, gridb2); + IdxRangeB1B2 gridb(gridb1, gridb2); - DFieldB1B2 results(gridb); + DFieldMemB1B2 results(gridb); quad_batched_operator( Kokkos::DefaultExecutionSpace(), - results.span_view(), + get_field(results), KOKKOS_LAMBDA(IdxXB1YB2 ibx) { - double const x = ddc::coordinate(ddc::select(ibx)); - double const y = ddc::coordinate(ddc::select(ibx)); - double const b1 = ddc::coordinate(ddc::select(ibx)); - double const b2 = ddc::coordinate(ddc::select(ibx)); + double const x = ddc::coordinate(ddc::select(ibx)); + double const y = ddc::coordinate(ddc::select(ibx)); + double const b1 = ddc::coordinate(ddc::select(ibx)); + double const b2 = ddc::coordinate(ddc::select(ibx)); return b1 * (x * y + 2 * x + b2 * y); }); auto results_host = ddc:: - create_mirror_view_and_copy(Kokkos::DefaultHostExecutionSpace(), results.span_view()); + create_mirror_view_and_copy(Kokkos::DefaultHostExecutionSpace(), get_field(results)); ddc::for_each(gridb, [&](IdxB1B2 ib) { - double const b1 = ddc::coordinate(ddc::select(ib)); - double const b2 = ddc::coordinate(ddc::select(ib)); + double const b1 = ddc::coordinate(ddc::select(ib)); + double const b2 = ddc::coordinate(ddc::select(ib)); double const exact = b1 * (y_max * y_max diff --git a/tests/quadrature/neumann_quadrature_spline.cpp b/tests/quadrature/neumann_quadrature_spline.cpp index 0fabc4ae5..b29ff2f4c 100644 --- a/tests/quadrature/neumann_quadrature_spline.cpp +++ b/tests/quadrature/neumann_quadrature_spline.cpp @@ -14,7 +14,7 @@ struct X static bool constexpr PERIODIC = false; }; -using CoordX = ddc::Coordinate; +using CoordX = Coord; struct BSplinesX : ddc::UniformBSplines { @@ -25,7 +25,7 @@ auto constexpr SplineXBoundary = ddc::BoundCond::HERMITE; using SplineInterpPointsX = ddc::KnotsAsInterpolationPoints; -struct IDimX : SplineInterpPointsX::interpolation_discrete_dimension_type +struct GridX : SplineInterpPointsX::interpolation_discrete_dimension_type { }; @@ -33,38 +33,39 @@ using SplineXBuilder_1d = ddc::SplineBuilder< Kokkos::DefaultHostExecutionSpace, Kokkos::DefaultHostExecutionSpace::memory_space, BSplinesX, - IDimX, + GridX, SplineXBoundary, SplineXBoundary, ddc::SplineSolver::LAPACK, - IDimX>; + GridX>; -using IVectX = ddc::DiscreteVector; -using IDomainX = ddc::DiscreteDomain; +using IdxStepX = IdxStep; +using IdxRangeX = IdxRange; -using DFieldX = ddc::Chunk; +using DFieldMemX = host_t>; TEST(NeumannSplineUniformQuadrature1D, ExactForConstantFunc) { CoordX const x_min(0.0); CoordX const x_max(M_PI); - IVectX const x_size(10); + IdxStepX const x_size(10); // Creating mesh & supports ddc::init_discrete_space(x_min, x_max, x_size); - ddc::init_discrete_space(SplineInterpPointsX::get_sampling()); - ddc::DiscreteDomain gridx(SplineInterpPointsX::get_domain()); + ddc::init_discrete_space(SplineInterpPointsX::get_sampling()); + IdxRangeX gridx(SplineInterpPointsX::get_domain()); SplineXBuilder_1d const builder_x(gridx); - DFieldX const quadrature_coeffs_host = neumann_spline_quadrature_coefficients(gridx, builder_x); + DFieldMemX const quadrature_coeffs_host + = neumann_spline_quadrature_coefficients(gridx, builder_x); auto quadrature_coeffs = ddc::create_mirror_and_copy( Kokkos::DefaultExecutionSpace(), - quadrature_coeffs_host.span_view()); + get_field(quadrature_coeffs_host)); Quadrature const integrate(quadrature_coeffs.span_cview()); - device_t values(gridx); + device_t values(gridx); ddc::parallel_fill(Kokkos::DefaultExecutionSpace(), values, 1.0); double integral = integrate(Kokkos::DefaultExecutionSpace(), values.span_cview()); double expected_val = x_max - x_min; @@ -84,7 +85,7 @@ struct ComputeErrorTraits }; using GrevillePointsY = ddc:: KnotsAsInterpolationPoints; - struct IDimY : GrevillePointsY::interpolation_discrete_dimension_type + struct GridY : GrevillePointsY::interpolation_discrete_dimension_type { }; }; @@ -95,44 +96,44 @@ double compute_error(int n_elems) using DimY = typename ComputeErrorTraits::Y; using BSplinesY = typename ComputeErrorTraits::BSplinesY; using GrevillePointsY = typename ComputeErrorTraits::GrevillePointsY; - using IDimY = typename ComputeErrorTraits::IDimY; + using GridY = typename ComputeErrorTraits::GridY; using SplineYBuilder = ddc::SplineBuilder< Kokkos::DefaultHostExecutionSpace, Kokkos::DefaultHostExecutionSpace::memory_space, BSplinesY, - IDimY, + GridY, ddc::BoundCond::HERMITE, ddc::BoundCond::HERMITE, ddc::SplineSolver::LAPACK, - IDimY>; - using IDomainY = ddc::DiscreteDomain; - using DFieldY = device_t>; - using DSpanY = device_t>; + GridY>; + using IdxRangeY = IdxRange; + using DFieldMemY = FieldMem; + using DFieldY = Field; - ddc::Coordinate const y_min(-1.0); - ddc::Coordinate const y_max(1.0); + Coord const y_min(-1.0); + Coord const y_max(1.0); ddc::init_discrete_space(y_min, y_max, n_elems); - ddc::init_discrete_space(GrevillePointsY::template get_sampling()); - IDomainY const gridy(GrevillePointsY::template get_domain()); + ddc::init_discrete_space(GrevillePointsY::template get_sampling()); + IdxRangeY const gridy(GrevillePointsY::template get_domain()); SplineYBuilder const builder_y(gridy); - host_t const quadrature_coeffs_host + host_t const quadrature_coeffs_host = neumann_spline_quadrature_coefficients(gridy, builder_y); auto quadrature_coeffs = ddc::create_mirror_and_copy( Kokkos::DefaultExecutionSpace(), - quadrature_coeffs_host.span_view()); + get_field(quadrature_coeffs_host)); Quadrature const integrate(quadrature_coeffs.span_cview()); - DFieldY values_alloc(gridy); - DSpanY values = values_alloc.span_view(); + DFieldMemY values_alloc(gridy); + DFieldY values = get_field(values_alloc); ddc::parallel_for_each( Kokkos::DefaultExecutionSpace(), gridy, - KOKKOS_LAMBDA(ddc::DiscreteElement const idx) { + KOKKOS_LAMBDA(Idx const idx) { double x = ddc::coordinate(idx); values(idx) = (x + 1) * (x + 1) * (x + 1) * (x - 1) * (x - 1); }); diff --git a/tests/quadrature/quadrature_1d.cpp b/tests/quadrature/quadrature_1d.cpp index d8eca4000..f3eb3533f 100644 --- a/tests/quadrature/quadrature_1d.cpp +++ b/tests/quadrature/quadrature_1d.cpp @@ -10,22 +10,22 @@ namespace { -struct RDimXPeriod +struct DimXPeriod { static bool constexpr PERIODIC = true; }; -struct IDimXPeriod : ddc::NonUniformPointSampling +struct GridXPeriod : NonUniformGridBase { }; -using IDomXPeriod = ddc::DiscreteDomain; -using CoordXPeriod = ddc::Coordinate; +using IDomXPeriod = IdxRange; +using CoordXPeriod = Coord; TEST(TrapezoidUniformPeriodicQuadrature1D, ExactForConstantFunc) { CoordXPeriod const x_min(0.0); CoordXPeriod const x_max(M_PI); - ddc::DiscreteVector const x_size(100); + IdxStep const x_size(100); // Creating mesh & supports std::vector point_sampling; @@ -36,20 +36,20 @@ TEST(TrapezoidUniformPeriodicQuadrature1D, ExactForConstantFunc) point_sampling.push_back(x_min + k * dx); } - ddc::init_discrete_space(point_sampling); - ddc::DiscreteElement lbound(0); - ddc::DiscreteVector npoints(x_size); - ddc::DiscreteDomain gridx(lbound, npoints); + ddc::init_discrete_space(point_sampling); + Idx lbound(0); + IdxStep npoints(x_size); + IdxRange gridx(lbound, npoints); - ddc::Chunk const quadrature_coeffs_host( + host_t> const quadrature_coeffs_host( trapezoid_quadrature_coefficients(gridx)); auto quadrature_coeffs = ddc::create_mirror_and_copy( Kokkos::DefaultExecutionSpace(), - quadrature_coeffs_host.span_view()); + get_field(quadrature_coeffs_host)); Quadrature const integrate(quadrature_coeffs.span_cview()); - device_t> values(gridx); - ddc::parallel_fill(Kokkos::DefaultExecutionSpace(), values.span_view(), 1.0); + FieldMem values(gridx); + ddc::parallel_fill(Kokkos::DefaultExecutionSpace(), get_field(values), 1.0); double integral = integrate(Kokkos::DefaultExecutionSpace(), values.span_cview()); double expected_val = x_max - x_min; EXPECT_LE(abs(integral - expected_val), 1e-9); @@ -59,7 +59,7 @@ TEST(SimpsonUniformPeriodicQuadrature1D, ExactForConstantFunc) { CoordXPeriod const x_min(0.0); CoordXPeriod const x_max(M_PI); - ddc::DiscreteVector const x_size(100); + IdxStep const x_size(100); // Creating mesh & support std::vector point_sampling; @@ -69,20 +69,20 @@ TEST(SimpsonUniformPeriodicQuadrature1D, ExactForConstantFunc) point_sampling.push_back(x_min + k * dx); } - ddc::init_discrete_space(point_sampling); - ddc::DiscreteElement lbound(0); - ddc::DiscreteVector npoints(x_size); - ddc::DiscreteDomain gridx(lbound, npoints); + ddc::init_discrete_space(point_sampling); + Idx lbound(0); + IdxStep npoints(x_size); + IdxRange gridx(lbound, npoints); - ddc::Chunk const quadrature_coeffs_host( + host_t> const quadrature_coeffs_host( simpson_quadrature_coefficients_1d(gridx)); auto quadrature_coeffs = ddc::create_mirror_and_copy( Kokkos::DefaultExecutionSpace(), - quadrature_coeffs_host.span_view()); + get_field(quadrature_coeffs_host)); Quadrature const integrate(quadrature_coeffs.span_cview()); - device_t> values(gridx); + FieldMem values(gridx); - ddc::parallel_fill(Kokkos::DefaultExecutionSpace(), values.span_view(), 1.0); + ddc::parallel_fill(Kokkos::DefaultExecutionSpace(), get_field(values), 1.0); double integral = integrate(Kokkos::DefaultExecutionSpace(), values.span_cview()); double expected_val = x_max - x_min; @@ -96,7 +96,7 @@ struct ComputeErrorTraits { static bool constexpr PERIODIC = false; }; - struct IDimY : ddc::NonUniformPointSampling + struct GridY : NonUniformGridBase { }; }; @@ -108,14 +108,14 @@ template double compute_error(int n_elems, Method meth) { using DimY = typename ComputeErrorTraits::Y; - using IDimY = typename ComputeErrorTraits::IDimY; - using IDomainY = ddc::DiscreteDomain; - using DFieldY = device_t>; - using DSpanY = device_t>; - - ddc::Coordinate const y_min(0.0); - ddc::Coordinate const y_max(M_PI); - std::vector> point_sampling; + using GridY = typename ComputeErrorTraits::GridY; + using IdxRangeY = IdxRange; + using DFieldMemY = FieldMem; + using DFieldY = Field; + + Coord const y_min(0.0); + Coord const y_max(M_PI); + std::vector> point_sampling; double dy = (y_max - y_min) / n_elems; // Create & intialize Uniform mesh @@ -123,12 +123,12 @@ double compute_error(int n_elems, Method meth) point_sampling.push_back(y_min + k * dy); } - ddc::init_discrete_space(point_sampling); - ddc::DiscreteElement lbound(0); - ddc::DiscreteVector npoints(n_elems); - ddc::DiscreteDomain gridy(lbound, npoints); + ddc::init_discrete_space(point_sampling); + Idx lbound(0); + IdxStep npoints(n_elems); + IdxRange gridy(lbound, npoints); - host_t quadrature_coeffs_host; + host_t quadrature_coeffs_host; switch (meth) { case Method::TRAPEZ: quadrature_coeffs_host = trapezoid_quadrature_coefficients(gridy); @@ -138,16 +138,16 @@ double compute_error(int n_elems, Method meth) auto quadrature_coeffs = ddc::create_mirror_and_copy( Kokkos::DefaultExecutionSpace(), - quadrature_coeffs_host.span_view()); + get_field(quadrature_coeffs_host)); Quadrature const integrate(quadrature_coeffs.span_cview()); - DFieldY values_alloc(gridy); - DSpanY values = values_alloc.span_view(); + DFieldMemY values_alloc(gridy); + DFieldY values = get_field(values_alloc); ddc::parallel_for_each( Kokkos::DefaultExecutionSpace(), gridy, - KOKKOS_LAMBDA(ddc::DiscreteElement const idx) { + KOKKOS_LAMBDA(Idx const idx) { values(idx) = Kokkos::sin(ddc::coordinate(idx)); }); double integral = integrate(Kokkos::DefaultExecutionSpace(), values.span_cview()); diff --git a/tests/quadrature/quadrature_2d.cpp b/tests/quadrature/quadrature_2d.cpp index 120f38143..cbf9ed5de 100644 --- a/tests/quadrature/quadrature_2d.cpp +++ b/tests/quadrature/quadrature_2d.cpp @@ -21,51 +21,51 @@ struct Y static constexpr bool PERIODIC = false; }; -using CoordX = ddc::Coordinate; -using CoordY = ddc::Coordinate; +using CoordX = Coord; +using CoordY = Coord; -struct IDimX : ddc::UniformPointSampling +struct GridX : UniformGridBase { }; -struct IDimY : ddc::UniformPointSampling +struct GridY : UniformGridBase { }; -using IdxXY = ddc::DiscreteElement; +using IdxXY = Idx; -using IVectX = ddc::DiscreteVector; -using IVectY = ddc::DiscreteVector; +using IdxStepX = IdxStep; +using IdxStepY = IdxStep; -using IDomainX = ddc::DiscreteDomain; -using IDomainY = ddc::DiscreteDomain; -using IDomainXY = ddc::DiscreteDomain; +using IdxRangeX = IdxRange; +using IdxRangeY = IdxRange; +using IdxRangeXY = IdxRange; -using DFieldXY = device_t>; +using DFieldMemXY = FieldMem; TEST(TrapezoidUniformNonPeriodicQuadrature2D, ExactForConstantFunc) { CoordX const x_min(0.0); CoordX const x_max(M_PI); - IVectX const x_size(10); + IdxStepX const x_size(10); CoordY const y_min(0.0); CoordY const y_max(20.0); - IVectY const y_size(10); + IdxStepY const y_size(10); // Creating mesh & supports - IDomainX const gridx = ddc::init_discrete_space(IDimX::init(x_min, x_max, x_size)); - IDomainY const gridy = ddc::init_discrete_space(IDimY::init(y_min, y_max, y_size)); + IdxRangeX const gridx = ddc::init_discrete_space(GridX::init(x_min, x_max, x_size)); + IdxRangeY const gridy = ddc::init_discrete_space(GridY::init(y_min, y_max, y_size)); - IDomainXY const gridxy(gridx, gridy); + IdxRangeXY const gridxy(gridx, gridy); - host_t const quadrature_coeffs_host = trapezoid_quadrature_coefficients(gridxy); + host_t const quadrature_coeffs_host = trapezoid_quadrature_coefficients(gridxy); auto quadrature_coeffs = ddc::create_mirror_and_copy( Kokkos::DefaultExecutionSpace(), - quadrature_coeffs_host.span_view()); - Quadrature const integrate(quadrature_coeffs); + get_field(quadrature_coeffs_host)); + Quadrature const integrate(quadrature_coeffs); - DFieldXY values(gridxy); + DFieldMemXY values(gridxy); ddc::parallel_fill(Kokkos::DefaultExecutionSpace(), values, 1.0); double integral = integrate(Kokkos::DefaultExecutionSpace(), values.span_cview()); @@ -77,27 +77,27 @@ void integrated_function_operator() { CoordX x_min(0.0); CoordX x_max(3.0); - IVectX x_ncells(4); + IdxStepX x_ncells(4); CoordY y_min(4.0); CoordY y_max(8.0); - IVectY y_ncells(16); + IdxStepY y_ncells(16); - IDomainX gridx = ddc::init_discrete_space(IDimX::init(x_min, x_max, x_ncells)); + IdxRangeX gridx = ddc::init_discrete_space(GridX::init(x_min, x_max, x_ncells)); - IDomainY gridy = ddc::init_discrete_space(IDimY::init(y_min, y_max, y_ncells)); - IDomainXY gridxy(gridx, gridy); + IdxRangeY gridy = ddc::init_discrete_space(GridY::init(y_min, y_max, y_ncells)); + IdxRangeXY gridxy(gridx, gridy); - host_t quad_coeffs_host_second = trapezoid_quadrature_coefficients(gridxy); + host_t quad_coeffs_host_second = trapezoid_quadrature_coefficients(gridxy); auto quad_coeffs_second = ddc::create_mirror_view_and_copy( Kokkos::DefaultExecutionSpace(), - quad_coeffs_host_second.span_view()); + get_field(quad_coeffs_host_second)); Quadrature func_operator(quad_coeffs_second.span_cview()); double const integral = func_operator( Kokkos::DefaultExecutionSpace(), KOKKOS_LAMBDA(IdxXY ixy) { - double y = ddc::coordinate(ddc::select(ixy)); - double x = ddc::coordinate(ddc::select(ixy)); + double y = ddc::coordinate(ddc::select(ixy)); + double x = ddc::coordinate(ddc::select(ixy)); return x * y + 2; }); EXPECT_DOUBLE_EQ(integral, 132.); @@ -119,10 +119,10 @@ struct ComputeErrorTraits { static bool constexpr PERIODIC = false; }; - struct IDimX : ddc::UniformPointSampling + struct GridX : UniformGridBase { }; - struct IDimY : ddc::UniformPointSampling + struct GridY : UniformGridBase { }; }; @@ -132,43 +132,43 @@ double compute_error(int n_elems) { using DimX = typename ComputeErrorTraits::X; using DimY = typename ComputeErrorTraits::Y; - using IDimX = typename ComputeErrorTraits::IDimX; - using IDimY = typename ComputeErrorTraits::IDimY; - using IVectX = ddc::DiscreteVector; - using IVectY = ddc::DiscreteVector; - using IDomainX = ddc::DiscreteDomain; - using IDomainY = ddc::DiscreteDomain; - using IDomainXY = ddc::DiscreteDomain; - using DFieldXY = device_t>; - using DSpanXY = device_t>; - - ddc::Coordinate const x_min(0.0); - ddc::Coordinate const x_max(M_PI); - IVectX x_size(n_elems); - - ddc::Coordinate const y_min(0.0); - ddc::Coordinate const y_max(M_PI); - IVectY y_size(n_elems); - - IDomainX const gridx - = ddc::init_discrete_space(IDimX::template init(x_min, x_max, x_size)); - IDomainY const gridy - = ddc::init_discrete_space(IDimY::template init(y_min, y_max, y_size)); - IDomainXY const gridxy(gridx, gridy); - - host_t const quadrature_coeffs_host = trapezoid_quadrature_coefficients(gridxy); + using GridX = typename ComputeErrorTraits::GridX; + using GridY = typename ComputeErrorTraits::GridY; + using IdxStepX = IdxStep; + using IdxStepY = IdxStep; + using IdxRangeX = IdxRange; + using IdxRangeY = IdxRange; + using IdxRangeXY = IdxRange; + using DFieldMemXY = FieldMem; + using DFieldXY = Field; + + Coord const x_min(0.0); + Coord const x_max(M_PI); + IdxStepX x_size(n_elems); + + Coord const y_min(0.0); + Coord const y_max(M_PI); + IdxStepY y_size(n_elems); + + IdxRangeX const gridx + = ddc::init_discrete_space(GridX::template init(x_min, x_max, x_size)); + IdxRangeY const gridy + = ddc::init_discrete_space(GridY::template init(y_min, y_max, y_size)); + IdxRangeXY const gridxy(gridx, gridy); + + host_t const quadrature_coeffs_host = trapezoid_quadrature_coefficients(gridxy); auto quadrature_coeffs = ddc::create_mirror_and_copy( Kokkos::DefaultExecutionSpace(), - quadrature_coeffs_host.span_view()); - Quadrature const integrate(quadrature_coeffs); + get_field(quadrature_coeffs_host)); + Quadrature const integrate(quadrature_coeffs); - DFieldXY values_alloc(gridxy); - DSpanXY values = values_alloc.span_view(); + DFieldMemXY values_alloc(gridxy); + DFieldXY values = get_field(values_alloc); ddc::parallel_for_each( Kokkos::DefaultExecutionSpace(), gridxy, - KOKKOS_LAMBDA(ddc::DiscreteElement const idx) { + KOKKOS_LAMBDA(Idx const idx) { double const y_cos = Kokkos::cos(ddc::get(ddc::coordinate(idx))); values(idx) = sin(ddc::get(ddc::coordinate(idx))) * y_cos * y_cos; }); diff --git a/tests/quadrature/quadrature_spline.cpp b/tests/quadrature/quadrature_spline.cpp index ffb9d1055..44eae4ea9 100644 --- a/tests/quadrature/quadrature_spline.cpp +++ b/tests/quadrature/quadrature_spline.cpp @@ -21,7 +21,7 @@ struct X #endif }; -using CoordX = ddc::Coordinate; +using CoordX = Coord; struct BSplinesX : ddc::UniformBSplines { @@ -32,45 +32,46 @@ auto constexpr SplineXBoundary = X::PERIODIC ? ddc::BoundCond::PERIODIC : ddc::B using SplineInterpPointsX = ddc::GrevilleInterpolationPoints; -struct IDimX : SplineInterpPointsX::interpolation_discrete_dimension_type +struct GridX : SplineInterpPointsX::interpolation_discrete_dimension_type { }; -using IVectX = ddc::DiscreteVector; -using IDomainX = ddc::DiscreteDomain; -using DFieldX = device_t>; +using IdxStepX = IdxStep; +using IdxRangeX = IdxRange; +using DFieldMemX = FieldMem; TEST(SplineUniformQuadrature, ExactForConstantFunc) { CoordX const x_min(0.0); CoordX const x_max(M_PI); - IVectX const x_size(10); + IdxStepX const x_size(10); using SplineXBuilder = ddc::SplineBuilder< Kokkos::DefaultHostExecutionSpace, Kokkos::DefaultHostExecutionSpace::memory_space, BSplinesX, - IDimX, + GridX, SplineXBoundary, SplineXBoundary, ddc::SplineSolver::LAPACK, - IDimX>; + GridX>; ddc::init_discrete_space(x_min, x_max, x_size); - ddc::init_discrete_space(SplineInterpPointsX::get_sampling()); - IDomainX gridx(SplineInterpPointsX::get_domain()); + ddc::init_discrete_space(SplineInterpPointsX::get_sampling()); + IdxRangeX gridx(SplineInterpPointsX::get_domain()); SplineXBuilder const builder_x(gridx); - host_t const quadrature_coeffs_host = spline_quadrature_coefficients(gridx, builder_x); + host_t const quadrature_coeffs_host + = spline_quadrature_coefficients(gridx, builder_x); auto quadrature_coeffs = ddc::create_mirror_and_copy( Kokkos::DefaultExecutionSpace(), - quadrature_coeffs_host.span_view()); + get_field(quadrature_coeffs_host)); Quadrature const integrate(quadrature_coeffs.span_cview()); - DFieldX values(gridx); + DFieldMemX values(gridx); ddc::parallel_fill(Kokkos::DefaultExecutionSpace(), values, 1.0); double integral = integrate(Kokkos::DefaultExecutionSpace(), values.span_cview()); @@ -93,7 +94,7 @@ struct ComputeErrorTraits BSplinesY, ddc::BoundCond::GREVILLE, ddc::BoundCond::GREVILLE>; - struct IDimY : GrevillePointsY::interpolation_discrete_dimension_type + struct GridY : GrevillePointsY::interpolation_discrete_dimension_type { }; }; @@ -104,44 +105,45 @@ double compute_error(int n_elems) using Y = typename ComputeErrorTraits::Y; using BSplinesY = typename ComputeErrorTraits::BSplinesY; using GrevillePointsY = typename ComputeErrorTraits::GrevillePointsY; - using IDimY = typename ComputeErrorTraits::IDimY; + using GridY = typename ComputeErrorTraits::GridY; auto constexpr SplineYBoundary = ddc::BoundCond::GREVILLE; using SplineYBuilder = ddc::SplineBuilder< Kokkos::DefaultHostExecutionSpace, Kokkos::DefaultHostExecutionSpace::memory_space, BSplinesY, - IDimY, + GridY, SplineYBoundary, SplineYBoundary, ddc::SplineSolver::LAPACK, - IDimY>; - using IDomainY = ddc::DiscreteDomain; - using DFieldY = device_t>; - using DSpanY = device_t>; + GridY>; + using IdxRangeY = IdxRange; + using DFieldMemY = FieldMem; + using DFieldY = Field; - ddc::Coordinate const y_min(0.0); - ddc::Coordinate const y_max(M_PI); + Coord const y_min(0.0); + Coord const y_max(M_PI); ddc::init_discrete_space(y_min, y_max, n_elems); - ddc::init_discrete_space(GrevillePointsY::template get_sampling()); - IDomainY const gridy(GrevillePointsY::template get_domain()); + ddc::init_discrete_space(GrevillePointsY::template get_sampling()); + IdxRangeY const gridy(GrevillePointsY::template get_domain()); SplineYBuilder const builder_y(gridy); - host_t const quadrature_coeffs_host = spline_quadrature_coefficients(gridy, builder_y); + host_t const quadrature_coeffs_host + = spline_quadrature_coefficients(gridy, builder_y); auto quadrature_coeffs = ddc::create_mirror_and_copy( Kokkos::DefaultExecutionSpace(), - quadrature_coeffs_host.span_view()); + get_field(quadrature_coeffs_host)); Quadrature const integrate(quadrature_coeffs.span_cview()); - DFieldY values_alloc(gridy); - DSpanY values = values_alloc.span_view(); + DFieldMemY values_alloc(gridy); + DFieldY values = get_field(values_alloc); ddc::parallel_for_each( Kokkos::DefaultExecutionSpace(), gridy, - KOKKOS_LAMBDA(ddc::DiscreteElement const idx) { + KOKKOS_LAMBDA(Idx const idx) { values(idx) = Kokkos::sin(ddc::coordinate(idx)); }); double integral = integrate(Kokkos::DefaultExecutionSpace(), values);