diff --git a/.clang-tidy b/.clang-tidy index 669ed9bfb..a5f1e0fda 100644 --- a/.clang-tidy +++ b/.clang-tidy @@ -2,5 +2,5 @@ # # SPDX-License-Identifier: MIT --- -Checks: '-*,bugprone-reserved-identifier,hicpp-avoid-c-arrays,modernize-use-nullptr' +Checks: '-*,bugprone-reserved-identifier,hicpp-avoid-c-arrays,modernize-use-nullptr,readability-avoid-const-params-in-decls' WarningsAsErrors: '*' diff --git a/examples/non_uniform_heat_equation.cpp b/examples/non_uniform_heat_equation.cpp index 1fcdf5218..9f5be0999 100644 --- a/examples/non_uniform_heat_equation.cpp +++ b/examples/non_uniform_heat_equation.cpp @@ -17,8 +17,8 @@ //! [vector_generator] std::vector generate_random_vector( int n, - int lower_bound, - int higher_bound) + double lower_bound, + double higher_bound) { std::random_device rd; std::mt19937 gen(rd()); diff --git a/include/ddc/kernels/fft.hpp b/include/ddc/kernels/fft.hpp index 44bac1c89..74a01083d 100644 --- a/include/ddc/kernels/fft.hpp +++ b/include/ddc/kernels/fft.hpp @@ -564,6 +564,7 @@ void core( } Kokkos::parallel_for( + "ddc_fft_normalization", Kokkos::RangePolicy( execSpace, 0, diff --git a/include/ddc/kernels/splines.hpp b/include/ddc/kernels/splines.hpp index 1cae877fd..27595fecd 100644 --- a/include/ddc/kernels/splines.hpp +++ b/include/ddc/kernels/splines.hpp @@ -20,6 +20,7 @@ #include "splines/spline_evaluator_2d.hpp" #include "splines/splines_linear_problem.hpp" #include "splines/splines_linear_problem_2x2_blocks.hpp" +#include "splines/splines_linear_problem_3x3_blocks.hpp" #include "splines/splines_linear_problem_band.hpp" #include "splines/splines_linear_problem_dense.hpp" #include "splines/splines_linear_problem_maker.hpp" diff --git a/include/ddc/kernels/splines/greville_interpolation_points.hpp b/include/ddc/kernels/splines/greville_interpolation_points.hpp index 1689297e7..47af761ce 100644 --- a/include/ddc/kernels/splines/greville_interpolation_points.hpp +++ b/include/ddc/kernels/splines/greville_interpolation_points.hpp @@ -127,6 +127,8 @@ class GrevilleInterpolationPoints * when uniform splines are used with an odd degree and with boundary conditions which * do not introduce additional interpolation points. * + * @tparam Sampling The discrete dimension supporting the Greville points. + * * @returns The mesh of uniform Greville points. */ template < @@ -143,6 +145,8 @@ class GrevilleInterpolationPoints /** * Get the NonUniformPointSampling defining the Greville points. * + * @tparam Sampling The discrete dimension supporting the Greville points. + * * @returns The mesh of non-uniform Greville points. */ template < @@ -262,6 +266,8 @@ class GrevilleInterpolationPoints /** * Get the domain which gives us access to all of the Greville points. * + * @tparam Sampling The discrete dimension supporting the Greville points. + * * @returns The domain of the Greville points. */ template diff --git a/include/ddc/kernels/splines/knots_as_interpolation_points.hpp b/include/ddc/kernels/splines/knots_as_interpolation_points.hpp index 24c34f132..23de02c19 100644 --- a/include/ddc/kernels/splines/knots_as_interpolation_points.hpp +++ b/include/ddc/kernels/splines/knots_as_interpolation_points.hpp @@ -24,6 +24,10 @@ namespace ddc { * In the case of strongly non-uniform splines this choice may result in a less * well conditioned problem, however most mathematical stability results are proven * with this choice of interpolation points. + * + * @tparam BSplines The type of the uniform or non-uniform spline basis whose knots are used as interpolation points. + * @tparam BcXmin The lower boundary condition. + * @tparam BcXmin The upper boundary condition. */ template class KnotsAsInterpolationPoints diff --git a/include/ddc/kernels/splines/null_extrapolation_rule.hpp b/include/ddc/kernels/splines/null_extrapolation_rule.hpp index baa5bc6ab..c26c3039d 100644 --- a/include/ddc/kernels/splines/null_extrapolation_rule.hpp +++ b/include/ddc/kernels/splines/null_extrapolation_rule.hpp @@ -8,8 +8,16 @@ namespace ddc { +/** + * @brief A functor describing a null extrapolation boundary value for 1D spline evaluator. + */ struct NullExtrapolationRule { + /** + * @brief Evaluates the spline at a coordinate outside of the domain. + * + * @return A double with the value of the function outside the domain (here, 0.). + */ template KOKKOS_FUNCTION double operator()(CoordType, ChunkSpan) const { diff --git a/include/ddc/kernels/splines/spline_builder.hpp b/include/ddc/kernels/splines/spline_builder.hpp index 9f0332e9c..697deb5d3 100644 --- a/include/ddc/kernels/splines/spline_builder.hpp +++ b/include/ddc/kernels/splines/spline_builder.hpp @@ -406,17 +406,13 @@ class SplineBuilder ddc::ChunkSpan spline, ddc::ChunkSpan vals, - std::optional> const derivs_xmin + std::optional< + ddc::ChunkSpan> + derivs_xmin = std::nullopt, - std::optional> const derivs_xmax + std::optional< + ddc::ChunkSpan> + derivs_xmax = std::nullopt) const; private: @@ -770,6 +766,7 @@ operator()( auto derivs_xmin_values = *derivs_xmin; auto const dx_proxy = m_dx; ddc::parallel_for_each( + "ddc_splines_hermite_compute_lower_coefficients", exec_space(), batch_domain(), KOKKOS_LAMBDA(typename batch_domain_type::discrete_element_type j) { @@ -787,6 +784,7 @@ operator()( auto const& interp_size_proxy = interpolation_domain().extents(); auto const& nbasis_proxy = ddc::discrete_space().nbasis(); ddc::parallel_for_each( + "ddc_splines_fill_rhs", exec_space(), batch_domain(), KOKKOS_LAMBDA(typename batch_domain_type::discrete_element_type j) { @@ -807,6 +805,7 @@ operator()( auto derivs_xmax_values = *derivs_xmax; auto const dx_proxy = m_dx; ddc::parallel_for_each( + "ddc_splines_hermite_compute_upper_coefficients", exec_space(), batch_domain(), KOKKOS_LAMBDA(typename batch_domain_type::discrete_element_type j) { @@ -820,12 +819,13 @@ operator()( } // TODO : Consider optimizing - // Allocate and fill a transposed version of spline in order to get dimension of interest as last dimension (optimal for GPU, necessary for Ginkgo) + // Allocate and fill a transposed version of spline in order to get dimension of interest as last dimension (optimal for GPU, necessary for Ginkgo). Also select only relevant rows in case of periodic boundaries ddc::Chunk spline_tr_alloc( batched_spline_tr_domain(), ddc::KokkosAllocator()); ddc::ChunkSpan spline_tr = spline_tr_alloc.span_view(); ddc::parallel_for_each( + "ddc_splines_transpose_rhs", exec_space(), batch_domain(), KOKKOS_LAMBDA(typename batch_domain_type::discrete_element_type const j) { @@ -841,8 +841,9 @@ operator()( batch_domain().size()); // Compute spline coef matrix->solve(bcoef_section); - // Transpose back spline_tr in spline + // Transpose back spline_tr into spline. ddc::parallel_for_each( + "ddc_splines_transpose_back_rhs", exec_space(), batch_domain(), KOKKOS_LAMBDA(typename batch_domain_type::discrete_element_type const j) { @@ -852,9 +853,10 @@ operator()( } }); - // Not sure yet of what this part do + // Duplicate the lower spline coefficients to the upper side in case of periodic boundaries if (bsplines_type::is_periodic()) { ddc::parallel_for_each( + "ddc_splines_periodic_rows_duplicate_rhs", exec_space(), batch_domain(), KOKKOS_LAMBDA(typename batch_domain_type::discrete_element_type const j) { diff --git a/include/ddc/kernels/splines/spline_builder_2d.hpp b/include/ddc/kernels/splines/spline_builder_2d.hpp index 480222f9f..187e9d7b7 100644 --- a/include/ddc/kernels/splines/spline_builder_2d.hpp +++ b/include/ddc/kernels/splines/spline_builder_2d.hpp @@ -350,53 +350,37 @@ class SplineBuilder2D ddc::ChunkSpan spline, ddc::ChunkSpan vals, - std::optional> const derivs_min1 + std::optional< + ddc::ChunkSpan> + derivs_min1 = std::nullopt, - std::optional> const derivs_max1 + std::optional< + ddc::ChunkSpan> + derivs_max1 = std::nullopt, - std::optional> const derivs_min2 + std::optional< + ddc::ChunkSpan> + derivs_min2 = std::nullopt, - std::optional> const derivs_max2 + std::optional< + ddc::ChunkSpan> + derivs_max2 = std::nullopt, - std::optional> const mixed_derivs_min1_min2 + std::optional< + ddc::ChunkSpan> + mixed_derivs_min1_min2 = std::nullopt, - std::optional> const mixed_derivs_max1_min2 + std::optional< + ddc::ChunkSpan> + mixed_derivs_max1_min2 = std::nullopt, - std::optional> const mixed_derivs_min1_max2 + std::optional< + ddc::ChunkSpan> + mixed_derivs_min1_max2 = std::nullopt, - std::optional> const mixed_derivs_max1_max2 + std::optional< + ddc::ChunkSpan> + mixed_derivs_max1_max2 = std::nullopt) const; }; diff --git a/include/ddc/kernels/splines/spline_evaluator.hpp b/include/ddc/kernels/splines/spline_evaluator.hpp index d85d64b1b..97458707a 100644 --- a/include/ddc/kernels/splines/spline_evaluator.hpp +++ b/include/ddc/kernels/splines/spline_evaluator.hpp @@ -15,10 +15,23 @@ namespace ddc { +/** + * @brief A class to evaluate, differentiate or integrate a spline function. + * + * A class which contains an operator () which can be used to evaluate, differentiate or integrate a spline function. + * + * @tparam ExecSpace The Kokkos execution space on which the spline evaluation is performed. + * @tparam MemorySpace The Kokkos memory space on which the data (spline coefficients and evaluation) is stored. + * @tparam BSplines The discrete dimension representing the B-splines. + * @tparam EvaluationMesh The discrete dimension on which evaluation points are defined. + * @tparam LeftExtrapolationRule The lower extrapolation rule type. + * @tparam RightExtrapolationRule The upper extrapolation rule type. + * @tparam IDimX A variadic template of all the discrete dimensions forming the full space (EvaluationMesh + batched dimensions). + */ template < class ExecSpace, class MemorySpace, - class BSplinesType, + class BSplines, class EvaluationMesh, class LeftExtrapolationRule, class RightExtrapolationRule, @@ -26,46 +39,69 @@ template < class SplineEvaluator { private: - // Tags to determine what to evaluate + /** + * @brief Tag to indicate that the value of the spline should be evaluated. + */ struct eval_type { }; + /** + * @brief Tag to indicate that derivative of the spline should be evaluated. + */ struct eval_deriv_type { }; - using tag_type = typename BSplinesType::tag_type; + using tag_type = typename BSplines::tag_type; public: + /// @brief The type of the Kokkos execution space used by this class. using exec_space = ExecSpace; + /// @brief The type of the Kokkos memory space used by this class. using memory_space = MemorySpace; - using bsplines_type = BSplinesType; - - using left_extrapolation_rule_type = LeftExtrapolationRule; - using right_extrapolation_rule_type = RightExtrapolationRule; - + /// @brief The type of the evaluation discrete dimension (discrete dimension of interest) used by this class. using evaluation_mesh_type = EvaluationMesh; + /// @brief The discrete dimension representing the B-splines. + using bsplines_type = BSplines; + + /// @brief The type of the domain for the 1D evaluation mesh used by this class. using evaluation_domain_type = ddc::DiscreteDomain; + /// @brief The type of the whole domain representing evaluation points. using batched_evaluation_domain_type = ddc::DiscreteDomain; + /// @brief The type of the 1D spline domain corresponding to the dimension of interest. using spline_domain_type = ddc::DiscreteDomain; + /** + * @brief The type of the batch domain (obtained by removing the dimension of interest + * from the whole domain). + */ using batch_domain_type = typename ddc::detail::convert_type_seq_to_discrete_domain, ddc::detail::TypeSeq>>; + /** + * @brief The type of the whole spline domain (cartesian product of 1D spline domain + * and batch domain) preserving the order of dimensions. + */ using batched_spline_domain_type = typename ddc::detail::convert_type_seq_to_discrete_domain, ddc::detail::TypeSeq, ddc::detail::TypeSeq>>; + /// @brief The type of the extrapolation rule at the lower boundary. + using left_extrapolation_rule_type = LeftExtrapolationRule; + + /// @brief The type of the extrapolation rule at the upper boundary. + using right_extrapolation_rule_type = RightExtrapolationRule; + private: LeftExtrapolationRule m_left_extrap_rule; @@ -105,6 +141,14 @@ class SplineEvaluator memory_space>>, "RightExtrapolationRule::operator() has to be callable with usual arguments."); + /** + * @brief Build a SplineEvaluator acting on batched_spline_domain. + * + * @param left_extrap_rule The extrapolation rule at the lower boundary. + * @param right_extrap_rule The extrapolation rule at the upper boundary. + * + * @see NullExtrapolationRule ConstantExtrapolationRule PeriodicExtrapolationRule + */ explicit SplineEvaluator( LeftExtrapolationRule const& left_extrap_rule, RightExtrapolationRule const& right_extrap_rule) @@ -113,26 +157,79 @@ class SplineEvaluator { } + /** + * @brief Copy-constructs. + * + * @param x A reference to another SplineEvaluator. + */ SplineEvaluator(SplineEvaluator const& x) = default; + /** + * @brief Move-constructs. + * + * @param x An rvalue to another SplineEvaluator. + */ SplineEvaluator(SplineEvaluator&& x) = default; + /// @brief Destructs ~SplineEvaluator() = default; + /** + * @brief Copy-assigns. + * + * @param x A reference to another SplineEvaluator. + * @return A reference to this object. + */ SplineEvaluator& operator=(SplineEvaluator const& x) = default; + /** + * @brief Move-assigns. + * + * @param x An rvalue to another SplineEvaluator. + * @return A reference to this object. + */ SplineEvaluator& operator=(SplineEvaluator&& x) = default; + /** + * @brief Get the lower extrapolation rule. + * + * Extrapolation rules are functors used to define the behavior of the SplineEvaluator out of the domain where the break points of the B-splines are defined. + * + * @return The lower extrapolation rule. + * + * @see NullExtrapolationRule ConstantExtrapolationRule PeriodicExtrapolationRule + */ left_extrapolation_rule_type left_extrapolation_rule() const { return m_left_extrap_rule; } + /** + * @brief Get the upper extrapolation rule. + * + * Extrapolation rules are functors used to define the behavior of the SplineEvaluator out of the domain where the break points of the B-splines are defined. + * + * @return The upper extrapolation rule. + * + * @see NullExtrapolationRule ConstantExtrapolationRule PeriodicExtrapolationRule + */ right_extrapolation_rule_type right_extrapolation_rule() const { return m_right_extrap_rule; } + /** + * @brief Evaluate 1D spline function (described by its spline coefficients) at a given coordinate. + * + * The spline coefficients represent a 1D spline function defined on a B-splines (basis splines). They can be obtained via various methods, such as using a SplineBuilder. + * + * Remark: calling SplineBuilder then SplineEvaluator corresponds to a spline interpolation. + * + * @param coord_eval The coordinate where the spline is evaluated. Note that only the component along the dimension of interest is used. + * @param spline_coef A ChunkSpan storing the 1D spline coefficients. + * + * @return The value of the spline function at the desired coordinate. + */ template KOKKOS_FUNCTION double operator()( ddc::Coordinate const& coord_eval, @@ -142,6 +239,26 @@ class SplineEvaluator return eval(coord_eval, spline_coef); } + /** + * @brief Evaluate spline function (described by its spline coefficients) on a mesh. + * + * The spline coefficients represent a spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder. + * + * This is not a multidimensional evaluation. This is a batched 1D evaluation. This means that for each slice of coordinates + * identified by a batch_domain_type::discrete_element_type, the evaluation is performed with the 1D set of + * spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * Remark: calling SplineBuilder then SplineEvaluator corresponds to a spline interpolation. + * + * @param[out] spline_eval The values of the spline function at the desired coordinates. For practical reasons those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. + * @param[in] coords_eval The coordinates where the spline is evaluated. Those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. Note that the coordinates of the + * points represented by this domain are unused and irrelevant (but the points themselves (DiscreteElement) are used to select + * the set of 1D spline coefficients retained to perform the evaluation). + * @param[in] spline_coef A ChunkSpan storing the spline coefficients. + */ template void operator()( ddc::ChunkSpan const @@ -158,6 +275,7 @@ class SplineEvaluator batch_domain_type const batch_domain(spline_eval.domain()); ddc::parallel_for_each( + "ddc_splines_evaluate", exec_space(), batch_domain, KOKKOS_CLASS_LAMBDA(typename batch_domain_type::discrete_element_type const j) { @@ -170,6 +288,17 @@ class SplineEvaluator }); } + /** + * @brief Differentiate 1D spline function (described by its spline coefficients) at a given coordinate. + * + * The spline coefficients represent a 1D spline function defined on a B-splines (basis splines). They can be + * obtained via various methods, such as using a SplineBuilder. + * + * @param coord_eval The coordinate where the spline is differentiated. Note that only the component along the dimension of interest is used. + * @param spline_coef A ChunkSpan storing the 1D spline coefficients. + * + * @return The derivative of the spline function at the desired coordinate. + */ template KOKKOS_FUNCTION double deriv( ddc::Coordinate const& coord_eval, @@ -179,6 +308,24 @@ class SplineEvaluator return eval_no_bc(coord_eval, spline_coef); } + /** + * @brief Differentiate spline function (described by its spline coefficients) on a mesh. + * + * The spline coefficients represent a spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder. + * + * The derivation is not performed in a multidimensional way (in any sense). This is a batched 1D derivation. + * This means that for each slice of coordinates identified by a batch_domain_type::discrete_element_type, + * the derivation is performed with the 1D set of spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * @param[out] spline_eval The derivatives of the spline function at the desired coordinates. For practical reasons those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. + * @param[in] coords_eval The coordinates where the spline is differentiated. Those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. Note that the coordinates of the + * points represented by this domain are unused and irrelevant (but the points themselves (DiscreteElement) are used to select + * the set of 1D spline coefficients retained to perform the evaluation). + * @param[in] spline_coef A ChunkSpan storing the spline coefficients. + */ template void deriv( ddc::ChunkSpan const @@ -195,6 +342,7 @@ class SplineEvaluator batch_domain_type const batch_domain(spline_eval.domain()); ddc::parallel_for_each( + "ddc_splines_differentiate", exec_space(), batch_domain, KOKKOS_CLASS_LAMBDA(typename batch_domain_type::discrete_element_type const j) { @@ -208,6 +356,19 @@ class SplineEvaluator }); } + /** @brief Perform batched 1D integrations of a spline function (described by its spline coefficients) along the dimension of interest and store results on a subdomain of batch_domain. + * + * The spline coefficients represent a spline function defined on a B-splines (basis splines). They can be obtained via the SplineBuilder. + * + * The integration is not performed in a multidimensional way (in any sense). This is a batched 1D integration. + * This means that for each element of integrals, the integration is performed with the 1D set of + * spline coefficients identified by the same DiscreteElement. + * + * @param[out] integrals The integrals of the spline function on the subdomain of batch_domain. For practical reasons those are + * stored in a ChunkSpan defined on a batch_domain_type. Note that the coordinates of the + * points represented by this domain are unused and irrelevant. + * @param[in] spline_coef A ChunkSpan storing the spline coefficients. + */ template void integrate( ddc::ChunkSpan const integrals, @@ -220,10 +381,12 @@ class SplineEvaluator ddc::KokkosAllocator()); ddc::ChunkSpan values = values_alloc.span_view(); Kokkos::parallel_for( + "ddc_splines_integrate_bsplines", Kokkos::RangePolicy(0, 1), KOKKOS_LAMBDA(int) { ddc::discrete_space().integrals(values); }); ddc::parallel_for_each( + "ddc_splines_integrate", exec_space(), batch_domain, KOKKOS_LAMBDA(typename batch_domain_type::discrete_element_type const j) { diff --git a/include/ddc/kernels/splines/spline_evaluator_2d.hpp b/include/ddc/kernels/splines/spline_evaluator_2d.hpp index fa77349a2..53c13b74c 100644 --- a/include/ddc/kernels/splines/spline_evaluator_2d.hpp +++ b/include/ddc/kernels/splines/spline_evaluator_2d.hpp @@ -16,15 +16,29 @@ namespace ddc { /** - * @brief Define an evaluator 2D on B-splines. + * @brief A class to evaluate, differentiate or integrate a 2D spline function. + * + * A class which contains an operator () which can be used to evaluate, differentiate or integrate a 2D spline function. + * + * @tparam ExecSpace The Kokkos execution space on which the spline evaluation is performed. + * @tparam MemorySpace The Kokkos memory space on which the data (spline coefficients and evaluation) is stored. + * @tparam BSplines1 The discrete dimension representing the B-splines along the first dimension of interest. + * @tparam BSplines2 The discrete dimension representing the B-splines along the second dimension of interest. + * @tparam EvaluationMesh1 The first discrete dimension on which evaluation points are defined. + * @tparam EvaluationMesh2 The second discrete dimension on which evaluation points are defined. + * @tparam LeftExtrapolationRule1 The lower extrapolation rule type along first dimension of interest. + * @tparam RightExtrapolationRule1 The upper extrapolation rule type along first dimension of interest. + * @tparam LeftExtrapolationRule2 The lower extrapolation rule type along second dimension of interest. + * @tparam RightExtrapolationRule2 The upper extrapolation rule type along second dimension of interest. + * @tparam IDimX A variadic template of all the discrete dimensions forming the full space (EvaluationMesh1 + EvaluationMesh2 + batched dimensions). */ template < class ExecSpace, class MemorySpace, - class BSplinesType1, - class BSplinesType2, - class evaluation_mesh_type1, - class evaluation_mesh_type2, + class BSplines1, + class BSplines2, + class EvaluationMesh1, + class EvaluationMesh2, class LeftExtrapolationRule1, class RightExtrapolationRule1, class LeftExtrapolationRule2, @@ -47,44 +61,80 @@ class SplineEvaluator2D { }; - using tag_type1 = typename BSplinesType1::tag_type; - using tag_type2 = typename BSplinesType2::tag_type; + using tag_type1 = typename BSplines1::tag_type; + using tag_type2 = typename BSplines2::tag_type; public: + /// @brief The type of the Kokkos execution space used by this class. using exec_space = ExecSpace; + /// @brief The type of the Kokkos memory space used by this class. using memory_space = MemorySpace; - using bsplines_type1 = BSplinesType1; - using bsplines_type2 = BSplinesType2; + /// @brief The type of the first discrete dimension of interest used by this class. + using evaluation_mesh_type1 = EvaluationMesh1; - using left_extrapolation_rule_1_type = LeftExtrapolationRule1; - using right_extrapolation_rule_1_type = RightExtrapolationRule1; - using left_extrapolation_rule_2_type = LeftExtrapolationRule2; - using right_extrapolation_rule_2_type = RightExtrapolationRule2; + /// @brief The type of the second discrete dimension of interest used by this class. + using evaluation_mesh_type2 = EvaluationMesh2; + /// @brief The discrete dimension representing the B-splines along first dimension. + using bsplines_type1 = BSplines1; + + /// @brief The discrete dimension representing the B-splines along second dimension. + using bsplines_type2 = BSplines2; + + /// @brief The type of the domain for the 1D evaluation mesh along first dimension used by this class. using evaluation_domain_type1 = ddc::DiscreteDomain; + + /// @brief The type of the domain for the 1D evaluation mesh along second dimension used by this class. using evaluation_domain_type2 = ddc::DiscreteDomain; + + /// @brief The type of the domain for the 2D evaluation mesh used by this class. using evaluation_domain_type = ddc::DiscreteDomain; + /// @brief The type of the whole domain representing evaluation points. using batched_evaluation_domain_type = ddc::DiscreteDomain; + /// @brief The type of the 1D spline domain corresponding to the first dimension of interest. using spline_domain_type1 = ddc::DiscreteDomain; + + /// @brief The type of the 1D spline domain corresponding to the second dimension of interest. using spline_domain_type2 = ddc::DiscreteDomain; + + /// @brief The type of the 2D spline domain corresponding to the dimensions of interest. using spline_domain_type = ddc::DiscreteDomain; + /** + * @brief The type of the batch domain (obtained by removing the dimensions of interest + * from the whole domain). + */ using batch_domain_type = typename ddc::detail::convert_type_seq_to_discrete_domain, ddc::detail::TypeSeq>>; + /** + * @brief The type of the whole spline domain (cartesian product of 2D spline domain + * and batch domain) preserving the underlying memory layout (order of dimensions). + */ using batched_spline_domain_type = typename ddc::detail::convert_type_seq_to_discrete_domain, ddc::detail::TypeSeq, ddc::detail::TypeSeq>>; + /// @brief The type of the extrapolation rule at the lower boundary along the first dimension. + using left_extrapolation_rule_1_type = LeftExtrapolationRule1; + + /// @brief The type of the extrapolation rule at the upper boundary along the first dimension. + using right_extrapolation_rule_1_type = RightExtrapolationRule1; + + /// @brief The type of the extrapolation rule at the lower boundary along the second dimension. + using left_extrapolation_rule_2_type = LeftExtrapolationRule2; + + /// @brief The type of the extrapolation rule at the upper boundary along the second dimension. + using right_extrapolation_rule_2_type = RightExtrapolationRule2; private: LeftExtrapolationRule1 m_left_extrap_rule_1; @@ -163,22 +213,14 @@ class SplineEvaluator2D "with usual arguments."); /** - * @brief Instantiate an evaluator operator. - * - * @param[in] left_extrap_rule1 - * A SplineBoundaryValue2D object giving the value on the "left side" of the domain - * in the first dimension. - * @param[in] right_extrap_rule1 - * A SplineBoundaryValue2D object giving the value on the "right side" of the domain - * in the first dimension. - * @param[in] left_extrap_rule2 - * A SplineBoundaryValue2D object giving the value on the "left side" of the domain - * in the second dimension. - * @param[in] right_extrap_rule2 - * A SplineBoundaryValue2D object giving the value on the "right side" of the domain - * in the second dimension. - * - * @see SplineBoundaryValue2D + * @brief Build a SplineEvaluator2D acting on batched_spline_domain. + * + * @param left_extrap_rule1 The extrapolation rule at the lower boundary along the first dimension. + * @param right_extrap_rule1 The extrapolation rule at the upper boundary along the first dimension. + * @param left_extrap_rule2 The extrapolation rule at the lower boundary along the second dimension. + * @param right_extrap_rule2 The extrapolation rule at the upper boundary along the second dimension. + * + * @see NullExtrapolationRule ConstantExtrapolationRule PeriodicExtrapolationRule */ explicit SplineEvaluator2D( LeftExtrapolationRule1 const& left_extrap_rule1, @@ -193,76 +235,105 @@ class SplineEvaluator2D } /** - * @brief Instantiate a SplineEvaluator2D from another - * SplineEvaluator2D (lvalue). + * @brief Copy-constructs. * - * @param[in] x - * SplineEvaluator2D evaluator used to instantiate the new one. + * @param x A reference to another SplineEvaluator. */ SplineEvaluator2D(SplineEvaluator2D const& x) = default; /** - * @brief Instantiate a SplineEvaluator2D from another temporary - * SplineEvaluator2D (rvalue). + * @brief Move-constructs. * - * @param[in] x - * SplineEvaluator2D evaluator used to instantiate the new one. + * @param x An rvalue to another SplineEvaluator. */ SplineEvaluator2D(SplineEvaluator2D&& x) = default; + /// @brief Destructs. ~SplineEvaluator2D() = default; /** - * @brief Assign a SplineEvaluator2D from another SplineEvaluator2D (lvalue). - * - * @param[in] x - * SplineEvaluator2D mapping used to assign. + * @brief Copy-assigns. * - * @return The SplineEvaluator2D assigned. + * @param x A reference to another SplineEvaluator. + * @return A reference to this object. */ SplineEvaluator2D& operator=(SplineEvaluator2D const& x) = default; /** - * @brief Assign a SplineEvaluator2D from another temporary SplineEvaluator2D (rvalue). + * @brief Move-assigns. * - * @param[in] x - * SplineEvaluator2D mapping used to assign. - * - * @return The SplineEvaluator2D assigned. + * @param x An rvalue to another SplineEvaluator. + * @return A reference to this object. */ SplineEvaluator2D& operator=(SplineEvaluator2D&& x) = default; - - + /** + * @brief Get the lower extrapolation rule along the first dimension. + * + * Extrapolation rules are functors used to define the behavior of the SplineEvaluator out of the domain where the break points of the B-splines are defined. + * + * @return The lower extrapolation rule along the first dimension. + * + * @see NullExtrapolationRule ConstantExtrapolationRule PeriodicExtrapolationRule + */ left_extrapolation_rule_1_type left_extrapolation_rule_dim_1() const { return m_left_extrap_rule_1; } + /** + * @brief Get the upper extrapolation rule along the first dimension. + * + * Extrapolation rules are functors used to define the behavior of the SplineEvaluator out of the domain where the break points of the B-splines are defined. + * + * @return The upper extrapolation rule along the first dimension. + * + * @see NullExtrapolationRule ConstantExtrapolationRule PeriodicExtrapolationRule + */ right_extrapolation_rule_1_type right_extrapolation_rule_dim_1() const { return m_right_extrap_rule_1; } + /** + * @brief Get the lower extrapolation rule along the second dimension. + * + * Extrapolation rules are functors used to define the behavior of the SplineEvaluator out of the domain where the break points of the B-splines are defined. + * + * @return The lower extrapolation rule along the second dimension. + * + * @see NullExtrapolationRule ConstantExtrapolationRule PeriodicExtrapolationRule + */ left_extrapolation_rule_2_type left_extrapolation_rule_dim_2() const { return m_left_extrap_rule_2; } + /** + * @brief Get the upper extrapolation rule along the second dimension. + * + * Extrapolation rules are functors used to define the behavior of the SplineEvaluator out of the domain where the break points of the B-splines are defined. + * + * @return The upper extrapolation rule along the second dimension. + * + * @see NullExtrapolationRule ConstantExtrapolationRule PeriodicExtrapolationRule + */ right_extrapolation_rule_2_type right_extrapolation_rule_dim_2() const { return m_right_extrap_rule_2; } /** - * @brief Get the value of the function on B-splines at the coordinate given. + * @brief Evaluate 2D spline function (described by its spline coefficients) at a given coordinate. * - * @param[in] coord_eval - * The 2D coordinate where we want to evaluate the function. - * @param[in] spline_coef - * The B-splines coefficients of the function we want to evaluate. + * The spline coefficients represent a 2D spline function defined on a B-splines (basis splines). They can be obtained via various methods, such as using a SplineBuilder2D. * - * @return A double containing the value of the function at the coordinate given. + * Remark: calling SplineBuilder2D then SplineEvaluator2D corresponds to a 2D spline interpolation. + * + * @param coord_eval The coordinate where the spline is evaluated. Note that only the components along the dimensions of interest are used. + * @param spline_coef A ChunkSpan storing the 2D spline coefficients. + * + * @return The value of the spline function at the desired coordinate. */ template KOKKOS_FUNCTION double operator()( @@ -273,6 +344,26 @@ class SplineEvaluator2D return eval(coord_eval, spline_coef); } + /** + * @brief Evaluate 2D spline function (described by its spline coefficients) on a mesh. + * + * The spline coefficients represent a 2D spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder2D. + * + * This is not a nD evaluation. This is a batched 2D evaluation. This means that for each slice of coordinates + * identified by a batch_domain_type::discrete_element_type, the evaluation is performed with the 2D set of + * spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * Remark: calling SplineBuilder2D then SplineEvaluator2D corresponds to a 2D spline interpolation. + * + * @param[out] spline_eval The values of the 2D spline function at the desired coordinates. For practical reasons those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. + * @param[in] coords_eval The coordinates where the spline is evaluated. Those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. Note that the coordinates of the + * points represented by this domain are unused and irrelevant (but the points themselves (DiscreteElement) are used to select + * the set of 2D spline coefficients retained to perform the evaluation). + * @param[in] spline_coef A ChunkSpan storing the 2D spline coefficients. + */ template void operator()( ddc::ChunkSpan const @@ -289,6 +380,7 @@ class SplineEvaluator2D evaluation_domain_type1 const evaluation_domain1(spline_eval.domain()); evaluation_domain_type2 const evaluation_domain2(spline_eval.domain()); ddc::parallel_for_each( + "ddc_splines_evaluate_2d", exec_space(), batch_domain, KOKKOS_CLASS_LAMBDA(typename batch_domain_type::discrete_element_type const j) { @@ -304,14 +396,15 @@ class SplineEvaluator2D } /** - * @brief Get the value of the derivative of the first dimension of the function on B-splines at the coordinate given. + * @brief Differentiate 2D spline function (described by its spline coefficients) at a given coordinate along first dimension of interest. * - * @param[in] coord_eval - * The 2D coordinate where we want to evaluate the derivative of the first dimension of the function. - * @param[in] spline_coef - * The B-splines coefficients of the function we want to evaluate. + * The spline coefficients represent a 2D spline function defined on a B-splines (basis splines). They can be + * obtained via various methods, such as using a SplineBuilder2D. * - * @return A double containing the value of the derivative of the first dimension of the function at the coordinate given. + * @param coord_eval The coordinate where the spline is differentiated. Note that only the components along the dimensions of interest are used. + * @param spline_coef A ChunkSpan storing the 2D spline coefficients. + * + * @return The derivative of the spline function at the desired coordinate. */ template KOKKOS_FUNCTION double deriv_dim_1( @@ -323,14 +416,15 @@ class SplineEvaluator2D } /** - * @brief Get the value of the derivative of the second dimension of the function on B-splines at the coordinate given. + * @brief Differentiate 2D spline function (described by its spline coefficients) at a given coordinate along second dimension of interest. * - * @param[in] coord_eval - * The 2D coordinate where we want to evaluate the derivative of the second dimension of the function. - * @param[in] spline_coef - * The B-splines coefficients of the function we want to evaluate. + * The spline coefficients represent a 2D spline function defined on a B-splines (basis splines). They can be + * obtained via various methods, such as using a SplineBuilder2D. + * + * @param coord_eval The coordinate where the spline is differentiated. Note that only the components along the dimensions of interest are used. + * @param spline_coef A ChunkSpan storing the 2D spline coefficients. * - * @return A double containing the value of the derivative of the second dimension of the function at the coordinate given. + * @return The derivative of the spline function at the desired coordinate. */ template KOKKOS_FUNCTION double deriv_dim_2( @@ -342,14 +436,15 @@ class SplineEvaluator2D } /** - * @brief Get the value of the cross derivative of the function on B-splines at the coordinate given. + * @brief Cross-differentiate 2D spline function (described by its spline coefficients) at a given coordinate. * - * @param[in] coord_eval - * The 2D coordinate where we want to evaluate the cross derivative of the function. - * @param[in] spline_coef - * The B-splines coefficients of the function we want to evaluate. + * The spline coefficients represent a 2D spline function defined on a B-splines (basis splines). They can be + * obtained via various methods, such as using a SplineBuilder2D. + * + * @param coord_eval The coordinate where the spline is differentiated. Note that only the components along the dimensions of interest are used. + * @param spline_coef A ChunkSpan storing the 2D spline coefficients. * - * @return A double containing the value of the cross derivative of the function at the coordinate given. + * @return The derivative of the spline function at the desired coordinate. */ template KOKKOS_FUNCTION double deriv_1_and_2( @@ -360,6 +455,19 @@ class SplineEvaluator2D return eval_no_bc(coord_eval, spline_coef); } + /** + * @brief Differentiate 2D spline function (described by its spline coefficients) at a given coordinate along a specified dimension of interest. + * + * The spline coefficients represent a 2D spline function defined on a B-splines (basis splines). They can be + * obtained via various methods, such as using a SplineBuilder2D. + * + * @tparam InterestDim Dimension along which differentiation is performed. + * + * @param coord_eval The coordinate where the spline is differentiated. Note that only the components along the dimensions of interest are used. + * @param spline_coef A ChunkSpan storing the 2D spline coefficients. + * + * @return The derivative of the spline function at the desired coordinate. + */ template KOKKOS_FUNCTION double deriv( ddc::Coordinate const& coord_eval, @@ -382,6 +490,22 @@ class SplineEvaluator2D } } + /** + * @brief Double-differentiate 2D spline function (described by its spline coefficients) at a given coordinate along specified dimensions of interest. + * + * The spline coefficients represent a 2D spline function defined on a B-splines (basis splines). They can be + * obtained via various methods, such as using a SplineBuilder2D. + * + * Note: double-differentiation other than cross-differentiation is not supported atm. See #440 + * + * @tparam InterestDim1 First dimension along which differentiation is performed. + * @tparam InterestDim2 Second dimension along which differentiation is performed. + * + * @param coord_eval The coordinate where the spline is double-differentiated. Note that only the components along the dimensions of interest are used. + * @param spline_coef A ChunkSpan storing the 2D spline coefficients. + * + * @return The derivative of the spline function at the desired coordinate. + */ template KOKKOS_FUNCTION double deriv2( ddc::Coordinate const& coord_eval, @@ -401,14 +525,24 @@ class SplineEvaluator2D } /** - * @brief Get the values of the derivative of the first dimension of the function on B-splines at the coordinates given. + * @brief Differentiate 2D spline function (described by its spline coefficients) on a mesh along first dimension of interest. * - * @param[out] spline_eval - * A ChunkSpan with the values of the derivative of the first dimension of the function at the coordinates given. - * @param[in] coords_eval - * A ChunkSpan with the 2D coordinates where we want to evaluate the derivative of the first dimension of the function. - * @param[in] spline_coef - * The B-splines coefficients of the function we want to evaluate. + * The spline coefficients represent a 2D spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder2D. + * + * This is not a nD evaluation. This is a batched 2D differentiation. + * This means that for each slice of coordinates identified by a batch_domain_type::discrete_element_type, + * the differentiation is performed with the 2D set of spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * @param[out] spline_eval The derivatives of the 2D spline function at the desired coordinates. For practical reasons those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. Note that the coordinates of the + * points represented by this domain are unused and irrelevant (but the points themselves (DiscreteElement) are used to select + * the set of 2D spline coefficients retained to perform the evaluation). + * @param[in] coords_eval The coordinates where the spline is differentiated. Those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. Note that the coordinates of the + * points represented by this domain are unused and irrelevant (but the points themselves (DiscreteElement) are used to select + * the set of 2D spline coefficients retained to perform the evaluation). + * @param[in] spline_coef A ChunkSpan storing the 2D spline coefficients. */ template void deriv_dim_1( @@ -426,6 +560,7 @@ class SplineEvaluator2D evaluation_domain_type1 const evaluation_domain1(spline_eval.domain()); evaluation_domain_type2 const evaluation_domain2(spline_eval.domain()); ddc::parallel_for_each( + "ddc_splines_differentiate_2d_dim_1", exec_space(), batch_domain, KOKKOS_CLASS_LAMBDA(typename batch_domain_type::discrete_element_type const j) { @@ -443,14 +578,22 @@ class SplineEvaluator2D } /** - * @brief Get the values of the derivative of the second dimension of the function on B-splines at the coordinates given. + * @brief Differentiate 2D spline function (described by its spline coefficients) on a mesh along second dimension of interest. * - * @param[out] spline_eval - * A ChunkSpan with the values of the derivative of the second dimension of the function at the coordinates given. - * @param[in] coords_eval - * A ChunkSpan with the 2D coordinates where we want to evaluate the derivative of the second dimension of the function. - * @param[in] spline_coef - * The B-splines coefficients of the function we want to evaluate. + * The spline coefficients represent a 2D spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder2D. + * + * This is not a nD differentiation. This is a batched 2D differentiation. + * This means that for each slice of coordinates identified by a batch_domain_type::discrete_element_type, + * the differentiation is performed with the 2D set of spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * @param[out] spline_eval The derivatives of the 2D spline function at the desired coordinates. For practical reasons those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. + * @param[in] coords_eval The coordinates where the spline is differentiated. Those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. Note that the coordinates of the + * points represented by this domain are unused and irrelevant (but the points themselves (DiscreteElement) are used to select + * the set of 2D spline coefficients retained to perform the evaluation). + * @param[in] spline_coef A ChunkSpan storing the 2D spline coefficients. */ template void deriv_dim_2( @@ -468,6 +611,7 @@ class SplineEvaluator2D evaluation_domain_type1 const evaluation_domain1(spline_eval.domain()); evaluation_domain_type2 const evaluation_domain2(spline_eval.domain()); ddc::parallel_for_each( + "ddc_splines_differentiate_2d_dim_2", exec_space(), batch_domain, KOKKOS_CLASS_LAMBDA(typename batch_domain_type::discrete_element_type const j) { @@ -485,14 +629,22 @@ class SplineEvaluator2D } /** - * @brief Get the values of the cross derivative of the function on B-splines at the coordinates given. + * @brief Cross-differentiate 2D spline function (described by its spline coefficients) on a mesh along dimensions of interest. * - * @param[out] spline_eval - * A ChunkSpan with the values of the cross derivative of the function at the coordinates given. - * @param[in] coords_eval - * A ChunkSpan with the 2D coordinates where we want to evaluate the cross derivative of the function. - * @param[in] spline_coef - * The B-splines coefficients of the function we want to evaluate. + * The spline coefficients represent a 2D spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder2D. + * + * This is not a nD cross-differentiation. This is a batched 2D cross-differentiation. + * This means that for each slice of coordinates identified by a batch_domain_type::discrete_element_type, + * the cross-differentiation is performed with the 2D set of spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * @param[out] spline_eval The cross-derivatives of the 2D spline function at the desired coordinates. For practical reasons those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. + * @param[in] coords_eval The coordinates where the spline is differentiated. Those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. Note that the coordinates of the + * points represented by this domain are unused and irrelevant (but the points themselves (DiscreteElement) are used to select + * the set of 2D spline coefficients retained to perform the evaluation). + * @param[in] spline_coef A ChunkSpan storing the 2D spline coefficients. */ template void deriv_1_and_2( @@ -510,6 +662,7 @@ class SplineEvaluator2D evaluation_domain_type1 const evaluation_domain1(spline_eval.domain()); evaluation_domain_type2 const evaluation_domain2(spline_eval.domain()); ddc::parallel_for_each( + "ddc_splines_cross_differentiate", exec_space(), batch_domain, KOKKOS_CLASS_LAMBDA(typename batch_domain_type::discrete_element_type const j) { @@ -526,6 +679,25 @@ class SplineEvaluator2D }); } + /** + * @brief Differentiate spline function (described by its spline coefficients) on a mesh along a specified dimension of interest. + * + * The spline coefficients represent a 2D spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder2D. + * + * This is not a nD evaluation. This is a batched 2D differentiation. + * This means that for each slice of coordinates identified by a batch_domain_type::discrete_element_type, + * the differentiation is performed with the 2D set of spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * @tparam InterestDim Dimension along which differentiation is performed. + * @param[out] spline_eval The derivatives of the 2D spline function at the desired coordinates. For practical reasons those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. + * @param[in] coords_eval The coordinates where the spline is differentiated. Those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. Note that the coordinates of the + * points represented by this domain are unused and irrelevant (but the points themselves (DiscreteElement) are used to select + * the set of 2D spline coefficients retained to perform the evaluation). + * @param[in] spline_coef A ChunkSpan storing the 2D spline coefficients. + */ template void deriv( ddc::ChunkSpan const @@ -554,6 +726,29 @@ class SplineEvaluator2D } } + /** + * @brief Double-differentiate 2D spline function (described by its spline coefficients) on a mesh along specified dimensions of interest. + * + * The spline coefficients represent a 2D spline function defined on a cartesian product of batch_domain and B-splines + * (basis splines). They can be obtained via various methods, such as using a SplineBuilder2D. + * + * This is not a nD evaluation. This is a batched 2D differentiation. + * This means that for each slice of coordinates identified by a batch_domain_type::discrete_element_type, + * the differentiation is performed with the 2D set of spline coefficients identified by the same batch_domain_type::discrete_element_type. + * + * Note: double-differentiation other than cross-differentiation is not supported atm. See #440 + * + * @tparam InterestDim1 First dimension along which differentiation is performed. + * @tparam InterestDim2 Second dimension along which differentiation is performed. + * + * @param[out] spline_eval The derivatives of the 2D spline function at the desired coordinates. For practical reasons those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. + * @param[in] coords_eval The coordinates where the spline is differentiated. Those are + * stored in a ChunkSpan defined on a batched_evaluation_domain_type. Note that the coordinates of the + * points represented by this domain are unused and irrelevant (but the points themselves (DiscreteElement) are used to select + * the set of 2D spline coefficients retained to perform the evaluation). + * @param[in] spline_coef A ChunkSpan storing the 2D spline coefficients. + */ template < class InterestDim1, class InterestDim2, @@ -584,12 +779,18 @@ class SplineEvaluator2D return deriv_1_and_2(spline_eval, coords_eval, spline_coef); } - /** - * @brief Get the the integral of the function on B-splines on the domain. - * @param[out] integrals - * The integrals of the function - * @param[in] spline_coef - * The B-splines coefficients of the function we want to integrate. + /** @brief Perform batched 2D integrations of a spline function (described by its spline coefficients) along the dimensions of interest and store results on a subdomain of batch_domain. + * + * The spline coefficients represent a 2D spline function defined on a B-splines (basis splines). They can be obtained via various methods, such as using a SplineBuilder2D. + * + * This is not a nD integration. This is a batched 2D integration. + * This means that for each element of integrals, the integration is performed with the 2D set of + * spline coefficients identified by the same DiscreteElement. + * + * @param[out] integrals The integrals of the 2D spline function on the subdomain of batch_domain. For practical reasons those are + * stored in a ChunkSpan defined on a batch_domain_type. Note that the coordinates of the + * points represented by this domain are unused and irrelevant. + * @param[in] spline_coef A ChunkSpan storing the 2D spline coefficients. */ template void integrate( @@ -607,6 +808,7 @@ class SplineEvaluator2D ddc::KokkosAllocator()); ddc::ChunkSpan values2 = values2_alloc.span_view(); Kokkos::parallel_for( + "ddc_splines_integrate_bsplines_2d", Kokkos::RangePolicy(0, 1), KOKKOS_LAMBDA(int) { ddc::discrete_space().integrals(values1); @@ -614,6 +816,7 @@ class SplineEvaluator2D }); ddc::parallel_for_each( + "ddc_splines_integrate_bsplines", exec_space(), batch_domain, KOKKOS_LAMBDA(typename batch_domain_type::discrete_element_type const j) { diff --git a/include/ddc/kernels/splines/splines_linear_problem_2x2_blocks.hpp b/include/ddc/kernels/splines/splines_linear_problem_2x2_blocks.hpp index c459166fb..fe605d71d 100644 --- a/include/ddc/kernels/splines/splines_linear_problem_2x2_blocks.hpp +++ b/include/ddc/kernels/splines/splines_linear_problem_2x2_blocks.hpp @@ -50,7 +50,7 @@ class SplinesLinearProblem2x2Blocks : public SplinesLinearProblem * @brief SplinesLinearProblem2x2Blocks constructor. * * @param mat_size The size of one of the dimensions of the square matrix. - * @param top_left_block A pointer toward the top-left SplinesLinearProblem. `setup_solver` must not have been called on `q`. + * @param top_left_block A pointer toward the top-left SplinesLinearProblem. `setup_solver` must not have been called on it. */ explicit SplinesLinearProblem2x2Blocks( std::size_t const mat_size, diff --git a/include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp b/include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp new file mode 100644 index 000000000..1dc93e99b --- /dev/null +++ b/include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp @@ -0,0 +1,185 @@ +// Copyright (C) The DDC development team, see COPYRIGHT.md file +// +// SPDX-License-Identifier: MIT + +#pragma once + +#include +#include +#include + +#include "splines_linear_problem.hpp" +#include "splines_linear_problem_2x2_blocks.hpp" + +namespace ddc::detail { + +/** + * @brief A 3x3-blocks linear problem dedicated to the computation of a spline approximation, + * with all blocks except center one being stored in dense format. + * + * A = | a | b | c | + * | d | Q | e | + * | f | g | h | + * + * The storage format is dense for all blocks except center one, whose storage format is determined by its type. + * + * The matrix itself and blocks a, Q and h are square (which fully determines the dimensions of the others). + * + * This class implements row & columns interchanges of the matrix and of multiple right-hand sides to restructure the + * 3x3-blocks linear problem into a 2x2-blocks linear problem, relying then on the operations implemented in SplinesLinearProblem2x2Blocks. + * + * @tparam ExecSpace The Kokkos::ExecutionSpace on which operations related to the matrix are supposed to be performed. + */ +template +class SplinesLinearProblem3x3Blocks : public SplinesLinearProblem2x2Blocks +{ +public: + using typename SplinesLinearProblem2x2Blocks::MultiRHS; + using SplinesLinearProblem2x2Blocks::size; + using SplinesLinearProblem2x2Blocks::solve; + using SplinesLinearProblem2x2Blocks::m_top_left_block; + +protected: + std::size_t m_top_size; + +public: + /** + * @brief SplinesLinearProblem3x3Blocks constructor. + * + * @param mat_size The size of one of the dimensions of the square matrix. + * @param top_size The size of one of the dimensions of the top-left square block. + * @param center_block A pointer toward the center SplinesLinearProblem. `setup_solver` must not have been called on it. + */ + explicit SplinesLinearProblem3x3Blocks( + std::size_t const mat_size, + std::size_t const top_size, + std::unique_ptr> center_block) + : SplinesLinearProblem2x2Blocks(mat_size, std::move(center_block)) + , m_top_size(top_size) + { + } + +private: + /// @brief Adjust indices, governs the row & columns interchanges to restructure the 3x3-blocks matrix into a 2x2-blocks matrix. + void adjust_indices(std::size_t& i, std::size_t& j) const + { + std::size_t const nq = m_top_left_block->size(); // size of the center block + + if (i < m_top_size) { + i += nq; + } else if (i < m_top_size + nq) { + i -= m_top_size; + } + + if (j < m_top_size) { + j += nq; + } else if (j < m_top_size + nq) { + j -= m_top_size; + } + } + +public: + double get_element(std::size_t i, std::size_t j) const override + { + adjust_indices(i, j); + return SplinesLinearProblem2x2Blocks::get_element(i, j); + } + + void set_element(std::size_t i, std::size_t j, double const aij) override + { + adjust_indices(i, j); + return SplinesLinearProblem2x2Blocks::set_element(i, j, aij); + } + +private: + /** + * @brief Perform row interchanges on multiple right-hand sides to get a 2-blocks structure (matching the requirements + * of the SplinesLinearProblem2x2Blocks solver). + * + * | b_top | | b_center | + * | b_center | -> | b_top | -- Considered as a + * | b_bottom | | b_bottom | -- single bottom block + * + * @param b The multiple right-hand sides. + */ + void interchange_rows_from_3_to_2_blocks_rhs(MultiRHS const b) const + { + std::size_t const nq = m_top_left_block->size(); // size of the center block + + // prevent Kokkos::deep_copy(b_top_dst, b_top) to be a deep_copy between overlapping allocations + assert(nq >= m_top_size); + + MultiRHS const b_top = Kokkos:: + subview(b, std::pair {0, m_top_size}, Kokkos::ALL); + MultiRHS const b_center = Kokkos:: + subview(b, + std::pair {m_top_size, m_top_size + nq}, + Kokkos::ALL); + + MultiRHS const b_center_dst + = Kokkos::subview(b, std::pair {0, nq}, Kokkos::ALL); + MultiRHS const b_top_dst = Kokkos:: + subview(b, std::pair {nq, nq + m_top_size}, Kokkos::ALL); + + MultiRHS const buffer = Kokkos::create_mirror(ExecSpace(), b_center); + + Kokkos::deep_copy(buffer, b_center); + Kokkos::deep_copy(b_top_dst, b_top); + Kokkos::deep_copy(b_center_dst, buffer); + } + + /** + * @brief Perform row interchanges on multiple right-hand sides to restore its 3-blocks structure. + * + * | b_center | | b_top | + * | b_top | -> | b_center | + * | b_bottom | | b_bottom | + * + * @param b The multiple right-hand sides. + */ + void interchange_rows_from_2_to_3_blocks_rhs(MultiRHS const b) const + { + std::size_t const nq = m_top_left_block->size(); // size of the center block + + // prevent Kokkos::deep_copy(b_top, b_top_src) to be a deep_copy between overlapping allocations + assert(nq >= m_top_size); + + MultiRHS const b_center_src + = Kokkos::subview(b, std::pair {0, nq}, Kokkos::ALL); + MultiRHS const b_top_src = Kokkos:: + subview(b, std::pair {nq, nq + m_top_size}, Kokkos::ALL); + + MultiRHS const b_top = Kokkos:: + subview(b, std::pair {0, m_top_size}, Kokkos::ALL); + MultiRHS const b_center = Kokkos:: + subview(b, + std::pair {m_top_size, m_top_size + nq}, + Kokkos::ALL); + + MultiRHS const buffer = Kokkos::create_mirror(ExecSpace(), b_center); + + Kokkos::deep_copy(buffer, b_center_src); + Kokkos::deep_copy(b_top, b_top_src); + Kokkos::deep_copy(b_center, buffer); + } + +public: + /** + * @brief Solve the multiple right-hand sides linear problem Ax=b or its transposed version A^tx=b inplace. + * + * Perform row interchanges on multiple right-hand sides to obtain a 2x2-blocks linear problem and call the SplinesLinearProblem2x2Blocks solver. + * + * @param[in, out] b A 2D Kokkos::View storing the multiple right-hand sides of the problem and receiving the corresponding solution. + * @param transpose Choose between the direct or transposed version of the linear problem. + */ + void solve(MultiRHS const b, bool const transpose) const override + { + assert(b.extent(0) == size()); + + interchange_rows_from_3_to_2_blocks_rhs(b); + SplinesLinearProblem2x2Blocks::solve(b, transpose); + interchange_rows_from_2_to_3_blocks_rhs(b); + } +}; + +} // namespace ddc::detail diff --git a/include/ddc/kernels/splines/splines_linear_problem_maker.hpp b/include/ddc/kernels/splines/splines_linear_problem_maker.hpp index 3823ac89c..8fe8183d5 100644 --- a/include/ddc/kernels/splines/splines_linear_problem_maker.hpp +++ b/include/ddc/kernels/splines/splines_linear_problem_maker.hpp @@ -8,6 +8,7 @@ #include #include "splines_linear_problem_2x2_blocks.hpp" +#include "splines_linear_problem_3x3_blocks.hpp" #include "splines_linear_problem_band.hpp" #include "splines_linear_problem_dense.hpp" #include "splines_linear_problem_pds_band.hpp" @@ -65,15 +66,16 @@ class SplinesLinearProblemMaker } /** - * @brief Construct a 2x2-blocks linear problem with band "main" block (the one called - * Q in SplinesLinearProblem2x2Blocks). + * @brief Construct a 2x2-blocks or 3x3-blocks linear problem with band "main" block (the one called + * Q in SplinesLinearProblem2x2Blocks and SplinesLinearProblem3x3Blocks). * * @tparam the Kokkos::ExecutionSpace on which matrix-related operation will be performed. * @param n The size of one of the dimensions of the whole square matrix. * @param kl The number of subdiagonals in the band block. * @param ku The number of superdiagonals in the band block. * @param pds A boolean indicating if the band block is positive-definite symetric or not. - * @param bottom_right_size The size of one of the dimensions of the bottom-right block. + * @param bottom_right_size The size of one of the dimensions of the bottom-right square block. + * @param top_left_size The size of one of the dimensions of the top-left square block. * * @return The SplinesLinearProblem instance. */ @@ -84,13 +86,53 @@ class SplinesLinearProblemMaker int const kl, int const ku, bool const pds, - int const bottom_size) + int const bottom_right_size, + int const top_left_size = 0) { - int const top_size = n - bottom_size; - std::unique_ptr> top_left_block - = make_new_band(top_size, kl, ku, pds); + int const main_size = n - top_left_size - bottom_right_size; + std::unique_ptr> main_block + = make_new_band(main_size, kl, ku, pds); + if (top_left_size == 0) { + return std::make_unique< + SplinesLinearProblem2x2Blocks>(n, std::move(main_block)); + } return std::make_unique< - SplinesLinearProblem2x2Blocks>(n, std::move(top_left_block)); + SplinesLinearProblem3x3Blocks>(n, top_left_size, std::move(main_block)); + } + + /** + * @brief Construct a 2x2-blocks linear problem with band "main" block (the one called + * Q in SplinesLinearProblem2x2Blocks) and other blocks containing the "periodic parts" of + * a periodic band matrix. + * + * It simply calls make_new_block_matrix_with_band_main_block with bottom_size being + * max(kl, ku) (except if the allocation would be higher than instantiating a SplinesLinearProblemDense). + * + * @tparam the Kokkos::ExecutionSpace on which matrix-related operation will be performed. + * @param n The size of one of the dimensions of the whole square matrix. + * @param kl The number of subdiagonals in the band block. + * @param ku The number of superdiagonals in the band block. + * @param pds A boolean indicating if the band block is positive-definite symetric or not. + * + * @return The SplinesLinearProblem instance. + */ + template + static std::unique_ptr> make_new_periodic_band_matrix( + int const n, + int const kl, + int const ku, + bool const pds) + { + assert(kl < n); + assert(ku < n); + int const bottom_size = std::max(kl, ku); + int const top_size = n - bottom_size; + + if (bottom_size * (n + top_size) + (2 * kl + ku + 1) * top_size >= n * n) { + return std::make_unique>(n); + } + + return make_new_block_matrix_with_band_main_block(n, kl, ku, pds, bottom_size); } /** diff --git a/tests/splines/matrix.cpp b/tests/splines/matrix.cpp index a1e8b1fec..19f5d6ae0 100644 --- a/tests/splines/matrix.cpp +++ b/tests/splines/matrix.cpp @@ -229,7 +229,7 @@ TEST(Matrix, 2x2Blocks) std::size_t const k = 10; std::unique_ptr> top_left_block = std::make_unique< - ddc::detail::SplinesLinearProblemDense>(3); + ddc::detail::SplinesLinearProblemDense>(7); std::unique_ptr> matrix = std::make_unique>(N, std::move(top_left_block)); @@ -248,6 +248,31 @@ TEST(Matrix, 2x2Blocks) solve_and_validate(*matrix); } +TEST(Matrix, 3x3Blocks) +{ + std::size_t const N = 10; + std::size_t const k = 10; + std::unique_ptr> center_block + = std::make_unique< + ddc::detail::SplinesLinearProblemDense>(N - 5); + std::unique_ptr> matrix + = std::make_unique>(N, 2, std::move(center_block)); + + // Build a non-symmetric full-rank matrix (without zero) + for (std::size_t i(0); i < N; ++i) { + std::cout << i; + matrix->set_element(i, i, 3. / 4 * ((N + 1) * i + 1)); + for (std::size_t j(std::max(0, int(i) - int(k))); j < i; ++j) { + matrix->set_element(i, j, -(1. / 4) / k * (N * i + j + 1)); + } + for (std::size_t j(i + 1); j < std::min(N, i + k + 1); ++j) { + matrix->set_element(i, j, -(1. / 4) / k * (N * i + j + 1)); + } + } + + solve_and_validate(*matrix); +} class MatrixSizesFixture : public testing::TestWithParam> { @@ -339,6 +364,56 @@ TEST_P(MatrixSizesFixture, 2x2Blocks) solve_and_validate(*matrix); } +TEST_P(MatrixSizesFixture, 3x3Blocks) +{ + auto const [N, k] = GetParam(); + std::unique_ptr> matrix + = ddc::detail::SplinesLinearProblemMaker::make_new_block_matrix_with_band_main_block< + Kokkos::DefaultExecutionSpace>(N, k, k, false, 3, 2); + + // Build a non-symmetric full-rank band matrix + for (std::size_t i(0); i < N; ++i) { + matrix->set_element(i, i, 3. / 4 * ((N + 1) * i + 1)); + for (std::size_t j(std::max(0, int(i) - int(k))); j < i; ++j) { + matrix->set_element(i, j, -(1. / 4) / k * (N * i + j + 1)); + } + for (std::size_t j(i + 1); j < std::min(N, i + k + 1); ++j) { + matrix->set_element(i, j, -(1. / 4) / k * (N * i + j + 1)); + } + } + + solve_and_validate(*matrix); +} + +TEST_P(MatrixSizesFixture, PeriodicBand) +{ + auto const [N, k] = GetParam(); + + // Build a full-rank periodic band matrix permuted in such a way the band is shifted + for (std::ptrdiff_t s(-k + k / 2 + 1); s < static_cast(k - k / 2); ++s) { + std::unique_ptr> matrix + = ddc::detail::SplinesLinearProblemMaker::make_new_periodic_band_matrix< + Kokkos::DefaultExecutionSpace>( + N, + static_cast(k - s), + k + s, + false); + for (std::size_t i(0); i < N; ++i) { + for (std::size_t j(0); j < N; ++j) { + std::ptrdiff_t diag = ddc::detail:: + modulo(static_cast(j - i), static_cast(N)); + if (diag == s || diag == N + s) { + matrix->set_element(i, j, 2.0 * k + 1); + } else if (diag <= s + k || diag >= N + s - k) { + matrix->set_element(i, j, -1.); + } + } + } + + solve_and_validate(*matrix); + } +} + TEST_P(MatrixSizesFixture, Sparse) { auto const [N, k] = GetParam();