From d121aea520645df5a9805e058502cd267a16445d Mon Sep 17 00:00:00 2001 From: Baptiste Legouix Date: Fri, 21 Jun 2024 18:21:43 +0200 Subject: [PATCH 01/23] Rename "interpolation" to "evaluation" in SplineEvaluator(2D) (#436) --------- Co-authored-by: EmilyBourne Co-authored-by: Thomas Padioleau --- .../ddc/kernels/splines/spline_evaluator.hpp | 70 +++++---- .../kernels/splines/spline_evaluator_2d.hpp | 140 +++++++++--------- 2 files changed, 101 insertions(+), 109 deletions(-) diff --git a/include/ddc/kernels/splines/spline_evaluator.hpp b/include/ddc/kernels/splines/spline_evaluator.hpp index c35ea4148..d85d64b1b 100644 --- a/include/ddc/kernels/splines/spline_evaluator.hpp +++ b/include/ddc/kernels/splines/spline_evaluator.hpp @@ -19,7 +19,7 @@ template < class ExecSpace, class MemorySpace, class BSplinesType, - class InterpolationMesh, + class EvaluationMesh, class LeftExtrapolationRule, class RightExtrapolationRule, class... IDimX> @@ -47,23 +47,23 @@ class SplineEvaluator using left_extrapolation_rule_type = LeftExtrapolationRule; using right_extrapolation_rule_type = RightExtrapolationRule; - using interpolation_mesh_type = InterpolationMesh; + using evaluation_mesh_type = EvaluationMesh; - using interpolation_domain_type = ddc::DiscreteDomain; + using evaluation_domain_type = ddc::DiscreteDomain; - using batched_interpolation_domain_type = ddc::DiscreteDomain; + using batched_evaluation_domain_type = ddc::DiscreteDomain; using spline_domain_type = ddc::DiscreteDomain; using batch_domain_type = typename ddc::detail::convert_type_seq_to_discrete_domain, - ddc::detail::TypeSeq>>; + ddc::detail::TypeSeq>>; using batched_spline_domain_type = typename ddc::detail::convert_type_seq_to_discrete_domain, - ddc::detail::TypeSeq, + ddc::detail::TypeSeq, ddc::detail::TypeSeq>>; @@ -144,17 +144,17 @@ class SplineEvaluator template void operator()( - ddc::ChunkSpan const + ddc::ChunkSpan const spline_eval, ddc::ChunkSpan< ddc::Coordinate const, - batched_interpolation_domain_type, + batched_evaluation_domain_type, Layout2, memory_space> const coords_eval, ddc::ChunkSpan const spline_coef) const { - interpolation_domain_type const interpolation_domain(spline_eval.domain()); + evaluation_domain_type const evaluation_domain(spline_eval.domain()); batch_domain_type const batch_domain(spline_eval.domain()); ddc::parallel_for_each( @@ -164,7 +164,7 @@ class SplineEvaluator const auto spline_eval_1D = spline_eval[j]; const auto coords_eval_1D = coords_eval[j]; const auto spline_coef_1D = spline_coef[j]; - for (auto const i : interpolation_domain) { + for (auto const i : evaluation_domain) { spline_eval_1D(i) = eval(coords_eval_1D(i), spline_coef_1D); } }); @@ -181,17 +181,17 @@ class SplineEvaluator template void deriv( - ddc::ChunkSpan const + ddc::ChunkSpan const spline_eval, ddc::ChunkSpan< ddc::Coordinate const, - batched_interpolation_domain_type, + batched_evaluation_domain_type, Layout2, memory_space> const coords_eval, ddc::ChunkSpan const spline_coef) const { - interpolation_domain_type const interpolation_domain(spline_eval.domain()); + evaluation_domain_type const evaluation_domain(spline_eval.domain()); batch_domain_type const batch_domain(spline_eval.domain()); ddc::parallel_for_each( @@ -201,7 +201,7 @@ class SplineEvaluator const auto spline_eval_1D = spline_eval[j]; const auto coords_eval_1D = coords_eval[j]; const auto spline_coef_1D = spline_coef[j]; - for (auto const i : interpolation_domain) { + for (auto const i : evaluation_domain) { spline_eval_1D(i) = eval_no_bc(coords_eval_1D(i), spline_coef_1D); } @@ -242,28 +242,27 @@ class SplineEvaluator ddc::ChunkSpan const spline_coef) const { - ddc::Coordinate - coord_eval_interpolation - = ddc::select( - coord_eval); + ddc::Coordinate + coord_eval_interest + = ddc::select(coord_eval); if constexpr (bsplines_type::is_periodic()) { - if (coord_eval_interpolation < ddc::discrete_space().rmin() - || coord_eval_interpolation > ddc::discrete_space().rmax()) { - coord_eval_interpolation -= Kokkos::floor( - (coord_eval_interpolation - - ddc::discrete_space().rmin()) - / ddc::discrete_space().length()) - * ddc::discrete_space().length(); + if (coord_eval_interest < ddc::discrete_space().rmin() + || coord_eval_interest > ddc::discrete_space().rmax()) { + coord_eval_interest -= Kokkos::floor( + (coord_eval_interest + - ddc::discrete_space().rmin()) + / ddc::discrete_space().length()) + * ddc::discrete_space().length(); } } else { - if (coord_eval_interpolation < ddc::discrete_space().rmin()) { - return m_left_extrap_rule(coord_eval_interpolation, spline_coef); + if (coord_eval_interest < ddc::discrete_space().rmin()) { + return m_left_extrap_rule(coord_eval_interest, spline_coef); } - if (coord_eval_interpolation > ddc::discrete_space().rmax()) { - return m_right_extrap_rule(coord_eval_interpolation, spline_coef); + if (coord_eval_interest > ddc::discrete_space().rmax()) { + return m_right_extrap_rule(coord_eval_interest, spline_coef); } } - return eval_no_bc(coord_eval_interpolation, spline_coef); + return eval_no_bc(coord_eval_interest, spline_coef); } template @@ -280,14 +279,13 @@ class SplineEvaluator double, std::experimental::extents> const vals(vals_ptr.data()); - ddc::Coordinate - coord_eval_interpolation - = ddc::select( - coord_eval); + ddc::Coordinate + coord_eval_interest + = ddc::select(coord_eval); if constexpr (std::is_same_v) { - jmin = ddc::discrete_space().eval_basis(vals, coord_eval_interpolation); + jmin = ddc::discrete_space().eval_basis(vals, coord_eval_interest); } else if constexpr (std::is_same_v) { - jmin = ddc::discrete_space().eval_deriv(vals, coord_eval_interpolation); + jmin = ddc::discrete_space().eval_deriv(vals, coord_eval_interest); } double y = 0.0; for (std::size_t i = 0; i < bsplines_type::degree() + 1; ++i) { diff --git a/include/ddc/kernels/splines/spline_evaluator_2d.hpp b/include/ddc/kernels/splines/spline_evaluator_2d.hpp index e33ced4c6..fa77349a2 100644 --- a/include/ddc/kernels/splines/spline_evaluator_2d.hpp +++ b/include/ddc/kernels/splines/spline_evaluator_2d.hpp @@ -23,8 +23,8 @@ template < class MemorySpace, class BSplinesType1, class BSplinesType2, - class interpolation_mesh_type1, - class interpolation_mesh_type2, + class evaluation_mesh_type1, + class evaluation_mesh_type2, class LeftExtrapolationRule1, class RightExtrapolationRule1, class LeftExtrapolationRule2, @@ -63,12 +63,12 @@ class SplineEvaluator2D using left_extrapolation_rule_2_type = LeftExtrapolationRule2; using right_extrapolation_rule_2_type = RightExtrapolationRule2; - using interpolation_domain_type1 = ddc::DiscreteDomain; - using interpolation_domain_type2 = ddc::DiscreteDomain; - using interpolation_domain_type - = ddc::DiscreteDomain; + using evaluation_domain_type1 = ddc::DiscreteDomain; + using evaluation_domain_type2 = ddc::DiscreteDomain; + using evaluation_domain_type + = ddc::DiscreteDomain; - using batched_interpolation_domain_type = ddc::DiscreteDomain; + using batched_evaluation_domain_type = ddc::DiscreteDomain; using spline_domain_type1 = ddc::DiscreteDomain; using spline_domain_type2 = ddc::DiscreteDomain; @@ -77,12 +77,12 @@ class SplineEvaluator2D using batch_domain_type = typename ddc::detail::convert_type_seq_to_discrete_domain, - ddc::detail::TypeSeq>>; + ddc::detail::TypeSeq>>; using batched_spline_domain_type = typename ddc::detail::convert_type_seq_to_discrete_domain, - ddc::detail::TypeSeq, + ddc::detail::TypeSeq, ddc::detail::TypeSeq>>; @@ -275,19 +275,19 @@ class SplineEvaluator2D template void operator()( - ddc::ChunkSpan const + ddc::ChunkSpan const spline_eval, ddc::ChunkSpan< ddc::Coordinate const, - batched_interpolation_domain_type, + batched_evaluation_domain_type, Layout2, memory_space> const coords_eval, ddc::ChunkSpan const spline_coef) const { batch_domain_type batch_domain(coords_eval.domain()); - interpolation_domain_type1 const interpolation_domain1(spline_eval.domain()); - interpolation_domain_type2 const interpolation_domain2(spline_eval.domain()); + evaluation_domain_type1 const evaluation_domain1(spline_eval.domain()); + evaluation_domain_type2 const evaluation_domain2(spline_eval.domain()); ddc::parallel_for_each( exec_space(), batch_domain, @@ -295,8 +295,8 @@ class SplineEvaluator2D const auto spline_eval_2D = spline_eval[j]; const auto coords_eval_2D = coords_eval[j]; const auto spline_coef_2D = spline_coef[j]; - for (auto const i1 : interpolation_domain1) { - for (auto const i2 : interpolation_domain2) { + for (auto const i1 : evaluation_domain1) { + for (auto const i2 : evaluation_domain2) { spline_eval_2D(i1, i2) = eval(coords_eval_2D(i1, i2), spline_coef_2D); } } @@ -369,16 +369,15 @@ class SplineEvaluator2D static_assert( std::is_same_v< InterestDim, - typename interpolation_mesh_type1:: - continuous_dimension_type> || std::is_same_v); + typename evaluation_mesh_type1:: + continuous_dimension_type> || std::is_same_v); if constexpr (std::is_same_v< InterestDim, - typename interpolation_mesh_type1::continuous_dimension_type>) { + typename evaluation_mesh_type1::continuous_dimension_type>) { return deriv_dim_1(coord_eval, spline_coef); } else if constexpr (std::is_same_v< InterestDim, - typename interpolation_mesh_type2:: - continuous_dimension_type>) { + typename evaluation_mesh_type2::continuous_dimension_type>) { return deriv_dim_2(coord_eval, spline_coef); } } @@ -392,12 +391,12 @@ class SplineEvaluator2D static_assert( (std::is_same_v< InterestDim1, - typename interpolation_mesh_type1:: - continuous_dimension_type> && std::is_same_v) + typename evaluation_mesh_type1:: + continuous_dimension_type> && std::is_same_v) || (std::is_same_v< InterestDim2, - typename interpolation_mesh_type1:: - continuous_dimension_type> && std::is_same_v)); + typename evaluation_mesh_type1:: + continuous_dimension_type> && std::is_same_v)); return deriv_1_and_2(coord_eval, spline_coef); } @@ -413,19 +412,19 @@ class SplineEvaluator2D */ template void deriv_dim_1( - ddc::ChunkSpan const + ddc::ChunkSpan const spline_eval, ddc::ChunkSpan< ddc::Coordinate const, - batched_interpolation_domain_type, + batched_evaluation_domain_type, Layout2, memory_space> const coords_eval, ddc::ChunkSpan const spline_coef) const { batch_domain_type batch_domain(coords_eval.domain()); - interpolation_domain_type1 const interpolation_domain1(spline_eval.domain()); - interpolation_domain_type2 const interpolation_domain2(spline_eval.domain()); + evaluation_domain_type1 const evaluation_domain1(spline_eval.domain()); + evaluation_domain_type2 const evaluation_domain2(spline_eval.domain()); ddc::parallel_for_each( exec_space(), batch_domain, @@ -433,8 +432,8 @@ class SplineEvaluator2D const auto spline_eval_2D = spline_eval[j]; const auto coords_eval_2D = coords_eval[j]; const auto spline_coef_2D = spline_coef[j]; - for (auto const i1 : interpolation_domain1) { - for (auto const i2 : interpolation_domain2) { + for (auto const i1 : evaluation_domain1) { + for (auto const i2 : evaluation_domain2) { spline_eval_2D(i1, i2) = eval_no_bc< eval_deriv_type, eval_type>(coords_eval_2D(i1, i2), spline_coef_2D); @@ -455,19 +454,19 @@ class SplineEvaluator2D */ template void deriv_dim_2( - ddc::ChunkSpan const + ddc::ChunkSpan const spline_eval, ddc::ChunkSpan< ddc::Coordinate const, - batched_interpolation_domain_type, + batched_evaluation_domain_type, Layout2, memory_space> const coords_eval, ddc::ChunkSpan const spline_coef) const { batch_domain_type batch_domain(coords_eval.domain()); - interpolation_domain_type1 const interpolation_domain1(spline_eval.domain()); - interpolation_domain_type2 const interpolation_domain2(spline_eval.domain()); + evaluation_domain_type1 const evaluation_domain1(spline_eval.domain()); + evaluation_domain_type2 const evaluation_domain2(spline_eval.domain()); ddc::parallel_for_each( exec_space(), batch_domain, @@ -475,8 +474,8 @@ class SplineEvaluator2D const auto spline_eval_2D = spline_eval[j]; const auto coords_eval_2D = coords_eval[j]; const auto spline_coef_2D = spline_coef[j]; - for (auto const i1 : interpolation_domain1) { - for (auto const i2 : interpolation_domain2) { + for (auto const i1 : evaluation_domain1) { + for (auto const i2 : evaluation_domain2) { spline_eval_2D(i1, i2) = eval_no_bc< eval_type, eval_deriv_type>(coords_eval_2D(i1, i2), spline_coef_2D); @@ -497,19 +496,19 @@ class SplineEvaluator2D */ template void deriv_1_and_2( - ddc::ChunkSpan const + ddc::ChunkSpan const spline_eval, ddc::ChunkSpan< ddc::Coordinate const, - batched_interpolation_domain_type, + batched_evaluation_domain_type, Layout2, memory_space> const coords_eval, ddc::ChunkSpan const spline_coef) const { batch_domain_type batch_domain(coords_eval.domain()); - interpolation_domain_type1 const interpolation_domain1(spline_eval.domain()); - interpolation_domain_type2 const interpolation_domain2(spline_eval.domain()); + evaluation_domain_type1 const evaluation_domain1(spline_eval.domain()); + evaluation_domain_type2 const evaluation_domain2(spline_eval.domain()); ddc::parallel_for_each( exec_space(), batch_domain, @@ -517,8 +516,8 @@ class SplineEvaluator2D const auto spline_eval_2D = spline_eval[j]; const auto coords_eval_2D = coords_eval[j]; const auto spline_coef_2D = spline_coef[j]; - for (auto const i1 : interpolation_domain1) { - for (auto const i2 : interpolation_domain2) { + for (auto const i1 : evaluation_domain1) { + for (auto const i2 : evaluation_domain2) { spline_eval_2D(i1, i2) = eval_no_bc< eval_deriv_type, eval_deriv_type>(coords_eval_2D(i1, i2), spline_coef_2D); @@ -529,11 +528,11 @@ class SplineEvaluator2D template void deriv( - ddc::ChunkSpan const + ddc::ChunkSpan const spline_eval, ddc::ChunkSpan< ddc::Coordinate const, - batched_interpolation_domain_type, + batched_evaluation_domain_type, Layout2, memory_space> const coords_eval, ddc::ChunkSpan const @@ -542,16 +541,15 @@ class SplineEvaluator2D static_assert( std::is_same_v< InterestDim, - typename interpolation_mesh_type1:: - continuous_dimension_type> || std::is_same_v); + typename evaluation_mesh_type1:: + continuous_dimension_type> || std::is_same_v); if constexpr (std::is_same_v< InterestDim, - typename interpolation_mesh_type1::continuous_dimension_type>) { + typename evaluation_mesh_type1::continuous_dimension_type>) { return deriv_dim_1(spline_eval, coords_eval, spline_coef); } else if constexpr (std::is_same_v< InterestDim, - typename interpolation_mesh_type2:: - continuous_dimension_type>) { + typename evaluation_mesh_type2::continuous_dimension_type>) { return deriv_dim_2(spline_eval, coords_eval, spline_coef); } } @@ -564,11 +562,11 @@ class SplineEvaluator2D class Layout3, class... CoordsDims> void deriv2( - ddc::ChunkSpan const + ddc::ChunkSpan const spline_eval, ddc::ChunkSpan< ddc::Coordinate const, - batched_interpolation_domain_type, + batched_evaluation_domain_type, Layout2, memory_space> const coords_eval, ddc::ChunkSpan const @@ -577,12 +575,12 @@ class SplineEvaluator2D static_assert( (std::is_same_v< InterestDim1, - typename interpolation_mesh_type1:: - continuous_dimension_type> && std::is_same_v) + typename evaluation_mesh_type1:: + continuous_dimension_type> && std::is_same_v) || (std::is_same_v< InterestDim2, - typename interpolation_mesh_type1:: - continuous_dimension_type> && std::is_same_v)); + typename evaluation_mesh_type1:: + continuous_dimension_type> && std::is_same_v)); return deriv_1_and_2(spline_eval, coords_eval, spline_coef); } @@ -656,8 +654,8 @@ class SplineEvaluator2D ddc::ChunkSpan const spline_coef) const { - using Dim1 = typename interpolation_mesh_type1::continuous_dimension_type; - using Dim2 = typename interpolation_mesh_type2::continuous_dimension_type; + using Dim1 = typename evaluation_mesh_type1::continuous_dimension_type; + using Dim2 = typename evaluation_mesh_type2::continuous_dimension_type; if constexpr (bsplines_type1::is_periodic()) { if (ddc::get(coord_eval) < ddc::discrete_space().rmin() || ddc::get(coord_eval) > ddc::discrete_space().rmax()) { @@ -698,8 +696,8 @@ class SplineEvaluator2D } return eval_no_bc( ddc::Coordinate< - typename interpolation_mesh_type1::continuous_dimension_type, - typename interpolation_mesh_type2::continuous_dimension_type>( + typename evaluation_mesh_type1::continuous_dimension_type, + typename evaluation_mesh_type2::continuous_dimension_type>( ddc::get(coord_eval), ddc::get(coord_eval)), spline_coef); @@ -742,28 +740,24 @@ class SplineEvaluator2D double, std::experimental::extents> const vals2(vals2_ptr.data()); - ddc::Coordinate - coord_eval_interpolation1 - = ddc::select( + ddc::Coordinate + coord_eval_interest1 + = ddc::select( coord_eval); - ddc::Coordinate - coord_eval_interpolation2 - = ddc::select( + ddc::Coordinate + coord_eval_interest2 + = ddc::select( coord_eval); if constexpr (std::is_same_v) { - jmin1 = ddc::discrete_space() - .eval_basis(vals1, coord_eval_interpolation1); + jmin1 = ddc::discrete_space().eval_basis(vals1, coord_eval_interest1); } else if constexpr (std::is_same_v) { - jmin1 = ddc::discrete_space() - .eval_deriv(vals1, coord_eval_interpolation1); + jmin1 = ddc::discrete_space().eval_deriv(vals1, coord_eval_interest1); } if constexpr (std::is_same_v) { - jmin2 = ddc::discrete_space() - .eval_basis(vals2, coord_eval_interpolation2); + jmin2 = ddc::discrete_space().eval_basis(vals2, coord_eval_interest2); } else if constexpr (std::is_same_v) { - jmin2 = ddc::discrete_space() - .eval_deriv(vals2, coord_eval_interpolation2); + jmin2 = ddc::discrete_space().eval_deriv(vals2, coord_eval_interest2); } double y = 0.0; From 9f9292d837f5b4d8ebee76b2f7c5e0e04eaa5176 Mon Sep 17 00:00:00 2001 From: Baptiste Legouix Date: Fri, 21 Jun 2024 18:47:03 +0200 Subject: [PATCH 02/23] Splines doc minor corrections (#495) * init * minor * Update include/ddc/kernels/splines/bsplines_uniform.hpp Co-authored-by: EmilyBourne --------- Co-authored-by: EmilyBourne --- .../kernels/splines/bsplines_non_uniform.hpp | 26 +++++++------- .../ddc/kernels/splines/bsplines_uniform.hpp | 34 +++++++++---------- .../ddc/kernels/splines/spline_builder.hpp | 10 +++--- .../ddc/kernels/splines/spline_builder_2d.hpp | 10 +++--- 4 files changed, 40 insertions(+), 40 deletions(-) diff --git a/include/ddc/kernels/splines/bsplines_non_uniform.hpp b/include/ddc/kernels/splines/bsplines_non_uniform.hpp index f54249aad..0c3cd4e84 100644 --- a/include/ddc/kernels/splines/bsplines_non_uniform.hpp +++ b/include/ddc/kernels/splines/bsplines_non_uniform.hpp @@ -155,9 +155,9 @@ class NonUniformBSplines : detail::NonUniformBSplinesBase template Impl(RandomIt breaks_begin, RandomIt breaks_end); - /** @brief Copy-constructs from another Impl with a different Kokkos memory space + /** @brief Copy-constructs from another Impl with a different Kokkos memory space. * - * @param impl A reference to the other Impl + * @param impl A reference to the other Impl. */ template explicit Impl(Impl const& impl) @@ -166,32 +166,32 @@ class NonUniformBSplines : detail::NonUniformBSplinesBase { } - /** @brief Copy-constructs + /** @brief Copy-constructs. * - * @param x A reference to another Impl + * @param x A reference to another Impl. */ Impl(Impl const& x) = default; - /** @brief Move-constructs + /** @brief Move-constructs. * - * @param x An rvalue to another Impl + * @param x An rvalue to another Impl. */ Impl(Impl&& x) = default; - /// @brief Destructs + /// @brief Destructs. ~Impl() = default; - /** @brief Copy-assigns + /** @brief Copy-assigns. * - * @param x A reference to another Impl - * @return A reference to the copied Impl + * @param x A reference to another Impl. + * @return A reference to the copied Impl. */ Impl& operator=(Impl const& x) = default; - /** @brief Move-assigns + /** @brief Move-assigns. * - * @param x An rvalue to another Impl - * @return A reference to the moved Impl + * @param x An rvalue to another Impl. + * @return A reference to this object. */ Impl& operator=(Impl&& x) = default; diff --git a/include/ddc/kernels/splines/bsplines_uniform.hpp b/include/ddc/kernels/splines/bsplines_uniform.hpp index f7da6b3b6..f15fc54f6 100644 --- a/include/ddc/kernels/splines/bsplines_uniform.hpp +++ b/include/ddc/kernels/splines/bsplines_uniform.hpp @@ -111,11 +111,11 @@ class UniformBSplines : detail::UniformBSplinesBase public: Impl() = default; - /** Constructs a spline basis (B-splines) with n equidistant knots over \f$[a, b]\f$ + /** Constructs a spline basis (B-splines) with n equidistant knots over \f$[a, b]\f$. * - * @param rmin the real ddc::coordinate of the first knot - * @param rmax the real ddc::coordinate of the last knot - * @param ncells the number of cells in the range [rmin, rmax] + * @param rmin The real ddc::coordinate of the first knot. + * @param rmax The real ddc::coordinate of the last knot. + * @param ncells The number of cells in the range [rmin, rmax]. */ explicit Impl(ddc::Coordinate rmin, ddc::Coordinate rmax, std::size_t ncells) { @@ -132,9 +132,9 @@ class UniformBSplines : detail::UniformBSplinesBase ddc::DiscreteVector(degree()))); } - /** @brief Copy-constructs from another Impl with a different Kokkos memory space + /** @brief Copy-constructs from another Impl with a different Kokkos memory space. * - * @param impl A reference to the other Impl + * @param impl A reference to the other Impl. */ template explicit Impl(Impl const& impl) @@ -143,32 +143,32 @@ class UniformBSplines : detail::UniformBSplinesBase { } - /** @brief Copy-constructs + /** @brief Copy-constructs. * - * @param x A reference to another Impl + * @param x A reference to another Impl. */ Impl(Impl const& x) = default; - /** @brief Move-constructs + /** @brief Move-constructs. * - * @param x An rvalue to another Impl + * @param x An rvalue to another Impl. */ Impl(Impl&& x) = default; - /// @brief Destructs + /// @brief Destructs. ~Impl() = default; - /** @brief Copy-assigns + /** @brief Copy-assigns. * - * @param x A reference to another Impl - * @return A reference to the copied Impl + * @param x A reference to another Impl. + * @return A reference to the copied Impl. */ Impl& operator=(Impl const& x) = default; - /** @brief Move-assigns + /** @brief Move-assigns. * - * @param x An rvalue to another Impl - * @return A reference to the moved Impl + * @param x An rvalue to another Impl. + * @return A reference to this object. */ Impl& operator=(Impl&& x) = default; diff --git a/include/ddc/kernels/splines/spline_builder.hpp b/include/ddc/kernels/splines/spline_builder.hpp index cf396aeaf..9f0332e9c 100644 --- a/include/ddc/kernels/splines/spline_builder.hpp +++ b/include/ddc/kernels/splines/spline_builder.hpp @@ -234,22 +234,22 @@ class SplineBuilder preconditionner_max_block_size); } - /// @brief Copy-constructor is deleted + /// @brief Copy-constructor is deleted. SplineBuilder(SplineBuilder const& x) = delete; - /** @brief Move-constructs + /** @brief Move-constructs. * * @param x An rvalue to another SplineBuilder. */ SplineBuilder(SplineBuilder&& x) = default; - /// @brief Destructs + /// @brief Destructs. ~SplineBuilder() = default; - /// @brief Copy-assignment is deleted + /// @brief Copy-assignment is deleted. SplineBuilder& operator=(SplineBuilder const& x) = delete; - /** @brief Move-assigns + /** @brief Move-assigns. * * @param x An rvalue to another SplineBuilder. * @return A reference to this object. diff --git a/include/ddc/kernels/splines/spline_builder_2d.hpp b/include/ddc/kernels/splines/spline_builder_2d.hpp index ebd0b8c6c..480222f9f 100644 --- a/include/ddc/kernels/splines/spline_builder_2d.hpp +++ b/include/ddc/kernels/splines/spline_builder_2d.hpp @@ -217,23 +217,23 @@ class SplineBuilder2D { } - /// @brief Copy-constructor is deleted + /// @brief Copy-constructor is deleted. SplineBuilder2D(SplineBuilder2D const& x) = delete; /** - * @brief Move-constructs + * @brief Move-constructs. * * @param x An rvalue to another SplineBuilder2D. */ SplineBuilder2D(SplineBuilder2D&& x) = default; - /// @brief Destructs + /// @brief Destructs. ~SplineBuilder2D() = default; - /// @brief Copy-assignment is deleted + /// @brief Copy-assignment is deleted. SplineBuilder2D& operator=(SplineBuilder2D const& x) = delete; - /** @brief Move-assigns + /** @brief Move-assigns. * * @param x An rvalue to another SplineBuilder. * @return A reference to this object. From 7e336b7e9c93e5bab1500f4f7994999b6bb0546d Mon Sep 17 00:00:00 2001 From: Thomas Padioleau Date: Mon, 24 Jun 2024 12:18:44 +0200 Subject: [PATCH 03/23] Update pages.yml --- .github/workflows/pages.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/pages.yml b/.github/workflows/pages.yml index bede61d6d..e88213ba6 100644 --- a/.github/workflows/pages.yml +++ b/.github/workflows/pages.yml @@ -161,6 +161,7 @@ jobs: mkdir docs_out chmod a+rwx docs_out docker run -v ${PWD}:/src ghcr.io/cexa-project/ddc/doxygen:${GITHUB_SHA:0:7} bash /src/run.sh + cat docs_out/doxygen.log - name: Publish site if: ${{ github.event_name == 'push' && github.ref_name == 'main' && needs.id_repo.outputs.in_base_repo == 'true' }} run: | From 47fbbf1ab6e64dcaf309417a8228a26ab1cf8eaf Mon Sep 17 00:00:00 2001 From: Baptiste Legouix Date: Mon, 24 Jun 2024 13:43:48 +0200 Subject: [PATCH 04/23] minor fix (#500) --- include/ddc/kernels/splines/splines_linear_problem_sparse.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/ddc/kernels/splines/splines_linear_problem_sparse.hpp b/include/ddc/kernels/splines/splines_linear_problem_sparse.hpp index 209e3d247..38ba8d28c 100644 --- a/include/ddc/kernels/splines/splines_linear_problem_sparse.hpp +++ b/include/ddc/kernels/splines/splines_linear_problem_sparse.hpp @@ -192,7 +192,7 @@ class SplinesLinearProblemSparse : public SplinesLinearProblem * * Removes the zeros from the CSR object and instantiate a Ginkgo solver. It also constructs a transposed version of the solver. * - * The stopping criterion is a reduction factor ||x||/||b||<1e-15 with 1000 maximum iterations. + * The stopping criterion is a reduction factor ||Ax-b||/||b||<1e-15 with 1000 maximum iterations. */ void setup_solver() override { From d1b51675bcfb4531fe2c218fa4917ea28d132e1b Mon Sep 17 00:00:00 2001 From: thierryantoun <152993060+thierryantoun@users.noreply.github.com> Date: Mon, 24 Jun 2024 14:12:59 +0200 Subject: [PATCH 05/23] DDC documentation for non uniform heat equation (#441) * Write the heat equation code with uniform space discretization * Format * Forgot to replace one [x][y] * Update examples/uniform_heat_equation.cpp Co-authored-by: Thomas Padioleau * Update examples/uniform_heat_equation.cpp Co-authored-by: Thomas Padioleau * Update examples/uniform_heat_equation.cpp Co-authored-by: Thomas Padioleau * Update examples/uniform_heat_equation.cpp Co-authored-by: Thomas Padioleau * Update examples/uniform_heat_equation.cpp Co-authored-by: Thomas Padioleau * Update examples/CMakeLists.txt Co-authored-by: Thomas Padioleau * Update examples/uniform_heat_equation.cpp Co-authored-by: Thomas Padioleau * Update examples/uniform_heat_equation.cpp Co-authored-by: Thomas Padioleau * Update examples/uniform_heat_equation.cpp Co-authored-by: Thomas Padioleau * Format * Adjust the commented example to the uniform discretisation case * Format * Forgot to add file * Adding constructors for non uniform domains and writting the non-uniform example with documentation * Correcting ghost points * Add documentation * Correction to doc file * Format document * Correction * Fix first_steps --------- Co-authored-by: Thomas Padioleau --- docs/CMakeLists.txt | 3 +- docs/first_steps.md | 37 +- docs/going_further.md | 123 +++++++ docs/heat_equation.md | 8 - docs/non_uniform_heat_equation.md | 8 + examples/CMakeLists.txt | 7 + examples/non_uniform_heat_equation.cpp | 401 +++++++++++++++++++++ examples/uniform_heat_equation.cpp | 8 + include/ddc/non_uniform_point_sampling.hpp | 78 ++++ 9 files changed, 663 insertions(+), 10 deletions(-) create mode 100644 docs/going_further.md delete mode 100644 docs/heat_equation.md create mode 100644 docs/non_uniform_heat_equation.md create mode 100644 examples/non_uniform_heat_equation.cpp diff --git a/docs/CMakeLists.txt b/docs/CMakeLists.txt index 76974a080..411e35490 100644 --- a/docs/CMakeLists.txt +++ b/docs/CMakeLists.txt @@ -59,7 +59,8 @@ set(DOXYGEN_IMAGE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/images") doxygen_add_docs(doc "${CMAKE_CURRENT_SOURCE_DIR}/About.md" "${CMAKE_CURRENT_SOURCE_DIR}/first_steps.md" + "${CMAKE_CURRENT_SOURCE_DIR}/going_further.md" "${CMAKE_CURRENT_SOURCE_DIR}/uniform_heat_equation.md" - "${CMAKE_CURRENT_SOURCE_DIR}/heat_equation.md" + "${CMAKE_CURRENT_SOURCE_DIR}/non_uniform_heat_equation.md" "${DDC_SOURCE_DIR}/include/ddc/" ALL) diff --git a/docs/first_steps.md b/docs/first_steps.md index 1b30af633..3d10b84b6 100644 --- a/docs/first_steps.md +++ b/docs/first_steps.md @@ -5,7 +5,7 @@ Copyright (C) The DDC development team, see COPYRIGHT.md file SPDX-License-Identifier: MIT --> -In \subpage heat_equation "examples/uniform_heat_equation.cpp" is a DDC example implementing a forward +In \subpage uniform_heat_equation "examples/uniform_heat_equation.cpp" is a DDC example implementing a forward finite-difference solver for the heat equation over a rectangle 2D domain with periodic boundary conditions. @@ -175,6 +175,33 @@ And we display the initial data. \snippet uniform_heat_equation.cpp initial-display + +\snippet uniform_heat_equation.cpp time iteration + +To display the data, a chunk is created on the host. + + +\snippet uniform_heat_equation.cpp host-chunk + +We deepcopy the data from the `ghosted_last_temp` chunk to `ghosted_temp` on the host. + +\snippet uniform_heat_equation.cpp boundary conditions + +\snippet uniform_heat_equation.cpp initial-deepcopy + +And we display the initial data. + +\snippet uniform_heat_equation.cpp initial-display + +For the numerical scheme, two chunkspans are created: ++ `next_temp` a span excluding ghosts of the temperature at the time-step we will build. ++ `last_temp` a read-only view of the temperature at the previous time-step.Note that *span_cview* returns a read-only ChunkSpan. + +\snippet uniform_heat_equation.cpp manipulated views + +We then solve the equation. + +\snippet uniform_heat_equation.cpp numerical scheme # Time loop \snippet uniform_heat_equation.cpp time iteration @@ -182,11 +209,16 @@ And we display the initial data. ## Periodic conditions +\snippet uniform_heat_equation.cpp output + \snippet uniform_heat_equation.cpp boundary conditions ## Numerical scheme + +\snippet uniform_heat_equation.cpp swap + For the numerical scheme, two chunkspans are created: + `next_temp` a span excluding ghosts of the temperature at the time-step we will build. + `last_temp` a read-only view of the temperature at the previous time-step.Note that *span_cview* returns a read-only ChunkSpan. @@ -195,4 +227,7 @@ For the numerical scheme, two chunkspans are created: We then solve the equation. + +\snippet uniform_heat_equation.cpp final output + \snippet uniform_heat_equation.cpp numerical scheme diff --git a/docs/going_further.md b/docs/going_further.md new file mode 100644 index 000000000..3feec26cf --- /dev/null +++ b/docs/going_further.md @@ -0,0 +1,123 @@ +# Commented example: the non uniform heat equation {#going_further} + + +In \subpage non_uniform_heat_equation "examples/non_uniform_heat_equation.cpp" is a DDC example implementing a forward +finite-difference solver for the heat equation over a rectangle 2D domain with periodic boundary +conditions and non-uniform space discretization. + +As usual, the file starts with a few includes that will be used in the code. + +\snippet non_uniform_heat_equation.cpp includes + +# Differences with the uniform problem resolution + +For non-uniform discretization, differences arise compared to uniform discretization, particularly in terms of dimension labeling, domain construction, and even in the resolution of the problem. + +## Dimensions naming + +Just like the uniform case, we start by defining types that we later use to name our +dimensions. + +\snippet non_uniform_heat_equation.cpp X-dimension + +However, we then define the discretization as non uniform. + +\snippet non_uniform_heat_equation.cpp X-discretization + +We do the same thing for the second dimension `Y`. + +\snippet non_uniform_heat_equation.cpp Y-space + +And for this case, we kept an uniform discretization for the time dimension. + +\snippet non_uniform_heat_equation.cpp time-space + +## Domains + +### Dimension X + +Just like the uniform case, we define the main characteristics of the domain. + +\snippet non_uniform_heat_equation.cpp main-start-x-parameters + +The first step to create a `DiscreteDomain` with non-uniform spatial discretization is to create a C++ iterator containing the actual coordinates of the points along each dimension. For tutorial purposes, we have constructed a function `generate_random_vector` that generates a vector with random data ranging from the lower bound to the higher one to populate our vectors. + +For the `X` dimension, we want to build 4 domains: +* `x_domain`: the main domain from the start (included) to end (included) but excluding "ghost" + points. +* `ghosted_x_domain`: the domain including all "ghost" points. +* `x_pre_ghost`: the "ghost" points that come before the main domain. +* `x_post_ghost`: the "ghost" points that come after the main domain. + +To do so, we have to create the iterator for the main domain. + +\snippet non_uniform_heat_equation.cpp iterator_main-domain + +For the initialization of ghost points in the non-uniform case, it was necessary to explicitly describe the position of the pre-ghost and post-ghost points. For the pre-ghost, it was necessary to start from the first coordinate along x and subtract the difference between the last and the penultimate coordinate of the x dimension to maintain periodicity. For the post-ghost point, take the last point and add the difference between the first and second points along the x dimension. + +\snippet non_uniform_heat_equation.cpp ghost_points_x + +Then we can create our 4 domains using the `init_discretization` +function with the `init_ghosted` function that takes the vectors as parameters. + +\snippet non_uniform_heat_equation.cpp build-domains + +**To summarize,** in this section we saw how to build a domain with non-uniform space discretization: ++ Create an iterator for your domain. ++ Call the `init_discretization` function with `init_ghosted` if you need ghost points or with `init` for a regular domain. + +### Dimension Y + +For the `Y` dimension, we do the same. We start by defining our 3 vectors. + +\snippet non_uniform_heat_equation.cpp Y-vectors + +And we use them to build the 4 domains in the same way as for the X dimension. + +\snippet non_uniform_heat_equation.cpp build-Y-domain + +### Time dimension + +Then we handle the domains for the simulated time dimension. We first give the simulated time at which to stard and end the simulation. + +\snippet non_uniform_heat_equation.cpp main-start-t-parameters + +The CFL conditions are more challenging to achieve for the case of non-uniform discretization. +Since the spatial steps are not uniform, we first need to find the maximum of the inverse of the square of the spatial step for the `X` and `Y` dimensions. And then we obtain the minimum value for `dt`. + +\snippet non_uniform_heat_equation.cpp CFL-condition + +We can calculate the number of time steps and build the `DiscreteDomain` for the time dimension. + +\snippet non_uniform_heat_equation.cpp time-domain + +## Time loop + +Allocation and initialization are the same as for the uniform case. Let's focus on resolving the numerical scheme. +The main difference in solving the numerical equation is that we need to account for the fact that the values of dx and dy on the left and right sides are different. We use the functions `distance_at_left` and `distance_at_right` to solve the equation. + +\snippet non_uniform_heat_equation.cpp numerical scheme + + + + + + + + + + + + + + + + + + + + diff --git a/docs/heat_equation.md b/docs/heat_equation.md deleted file mode 100644 index 555ac22d5..000000000 --- a/docs/heat_equation.md +++ /dev/null @@ -1,8 +0,0 @@ -# examples/heat_equation.cpp {#heat_equation} - - -\include{lineno} heat_equation.cpp diff --git a/docs/non_uniform_heat_equation.md b/docs/non_uniform_heat_equation.md new file mode 100644 index 000000000..9e2bf9a35 --- /dev/null +++ b/docs/non_uniform_heat_equation.md @@ -0,0 +1,8 @@ +# examples/non_uniform_heat_equation.cpp {#non_uniform_heat_equation} + + +\include{lineno} non_uniform_heat_equation.cpp \ No newline at end of file diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index c0392b47d..c7b19a0d6 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -6,11 +6,18 @@ add_executable(heat_equation heat_equation.cpp) target_link_libraries(heat_equation PUBLIC DDC::DDC) add_executable(uniform_heat_equation uniform_heat_equation.cpp) +add_executable(non_uniform_heat_equation non_uniform_heat_equation.cpp) + +target_link_libraries(uniform_heat_equation PUBLIC DDC::DDC) +target_link_libraries(non_uniform_heat_equation PUBLIC DDC::DDC) target_link_libraries(uniform_heat_equation PUBLIC DDC::DDC) + if("${DDC_BUILD_PDI_WRAPPER}") target_link_libraries(heat_equation PUBLIC DDC::PDI_Wrapper) target_link_libraries(uniform_heat_equation PUBLIC DDC::PDI_Wrapper) + target_link_libraries(non_uniform_heat_equation PUBLIC DDC::PDI_Wrapper) + endif() if("${DDC_BUILD_KERNELS_FFT}") add_executable(heat_equation_spectral heat_equation_spectral.cpp) diff --git a/examples/non_uniform_heat_equation.cpp b/examples/non_uniform_heat_equation.cpp new file mode 100644 index 000000000..1fcdf5218 --- /dev/null +++ b/examples/non_uniform_heat_equation.cpp @@ -0,0 +1,401 @@ +// Copyright (C) The DDC development team, see COPYRIGHT.md file +// +// SPDX-License-Identifier: MIT + +//! [includes] +#include +#include +#include +#include +#include +#include +#include + +#include +//! [includes] + +//! [vector_generator] +std::vector generate_random_vector( + int n, + int lower_bound, + int higher_bound) +{ + std::random_device rd; + std::mt19937 gen(rd()); + std::uniform_real_distribution + dis(lower_bound, higher_bound); + + std::vector vec(n); + vec[0] = lower_bound; + vec[n - 1] = higher_bound; + + for (int i = 1; i < vec.size() - 1; ++i) { + vec[i] = dis(gen); + } + + std::sort(vec.begin(), vec.end()); + return vec; +} +//! [vector_generator] + +//! [X-dimension] +struct X; +//! [X-dimension] + +//! [X-discretization] +struct DDimX : ddc::NonUniformPointSampling +{ +}; +//! [X-discretization] + +//! [Y-space] +struct Y; +struct DDimY : ddc::NonUniformPointSampling +{ +}; +//! [Y-space] + +//! [time-space] +struct T; +struct DDimT : ddc::UniformPointSampling +{ +}; +//! [time-space] + +//! [display] +/** A function to pretty print the temperature + * @tparam ChunkType The type of chunk span. This way the template parameters are avoided, + * should be deduced by the compiler. + * @param time The time at which the output is made. + * @param temp The temperature at this time-step. + */ +template +void display(double time, ChunkType temp) +{ + double const mean_temp = ddc::transform_reduce( + temp.domain(), + 0., + ddc::reducer::sum(), + temp) + / temp.domain().size(); + std::cout << std::fixed << std::setprecision(3); + std::cout << "At t = " << time << ",\n"; + std::cout << " * mean temperature = " << mean_temp << "\n"; + ddc::ChunkSpan temp_slice + = temp[ddc::get_domain(temp).front() + + ddc::get_domain(temp).size() / 2]; + std::cout << " * temperature[y:" + << ddc::get_domain(temp).size() / 2 << "] = {"; + ddc::for_each( + ddc::get_domain(temp), + [=](ddc::DiscreteElement const ix) { + std::cout << std::setw(6) << temp_slice(ix); + }); + std::cout << " }" << std::endl; + +#ifdef DDC_BUILD_PDI_WRAPPER + ddc::PdiEvent("display") + .with("temp", temp) + .and_with("mean_temp", mean_temp) + .and_with("temp_slice", temp_slice); +#endif +} +//! [display] + + +//! [main-start] +//! [main-start-x-parameters] +int main(int argc, char** argv) +{ +#ifdef DDC_BUILD_PDI_WRAPPER + auto pdi_conf = PC_parse_string(""); + PDI_init(pdi_conf); +#endif + Kokkos::ScopeGuard const kokkos_scope(argc, argv); + ddc::ScopeGuard const ddc_scope(argc, argv); + + + //! [parameters] + double const x_start = -1.; + double const x_end = 1.; + std::size_t const nb_x_points = 10; + double const kx = .01; + //! [main-start-x-parameters] + //! [main-start-y-parameters] + double const y_start = -1.; + double const y_end = 1.; + std::size_t const nb_y_points = 100; + double const ky = .002; + //! [main-start-y-parameters] + //! [main-start-t-parameters] + double const start_time = 0.; + double const end_time = 10.; + //! [main-start-t-parameters] + std::ptrdiff_t const t_output_period = 10; + //! [parameters] + + //! [iterator_main-domain] + std::vector x_domain_vect + = generate_random_vector(nb_x_points, x_start, x_end); + //! [iterator_main-domain] + + std::size_t size_x = x_domain_vect.size(); + + //! [ghost_points_x] + std::vector x_pre_ghost_vect { + x_domain_vect.front() + - (x_domain_vect.back() - x_domain_vect[size_x - 2])}; + + std::vector x_post_ghost_vect { + x_domain_vect.back() + + (x_domain_vect[1] - x_domain_vect.front())}; + //! [ghost_points_x] + + //! [build-domains] + auto const [x_domain, ghosted_x_domain, x_pre_ghost, x_post_ghost] + = ddc::init_discrete_space(DDimX::init_ghosted( + x_domain_vect, + x_pre_ghost_vect, + x_post_ghost_vect)); + //! [build-domains] + + ddc::DiscreteDomain const + x_domain_begin(x_domain.front(), x_post_ghost.extents()); + ddc::DiscreteDomain const x_domain_end( + x_domain.back() - x_pre_ghost.extents() + 1, + x_pre_ghost.extents()); + + //! [Y-vectors] + std::vector y_domain_vect + = generate_random_vector(nb_y_points, y_start, y_end); + + std::size_t size_y = y_domain_vect.size(); + + //! [ghost_points_y] + std::vector y_pre_ghost_vect { + y_domain_vect.front() + - (y_domain_vect.back() - y_domain_vect[size_y - 2])}; + std::vector y_post_ghost_vect { + y_domain_vect.back() + + (y_domain_vect[1] - y_domain_vect.front())}; + //! [ghost_points_y] + + //! [Y-vectors] + + //! [build-Y-domain] + auto const [y_domain, ghosted_y_domain, y_pre_ghost, y_post_ghost] + = ddc::init_discrete_space(DDimY::init_ghosted( + y_domain_vect, + y_pre_ghost_vect, + y_post_ghost_vect)); + //! [build-Y-domain] + + ddc::DiscreteDomain const + y_domain_begin(y_domain.front(), y_post_ghost.extents()); + + ddc::DiscreteDomain const y_domain_end( + y_domain.back() - y_pre_ghost.extents() + 1, + y_pre_ghost.extents()); + + //! [CFL-condition] + + double const invdx2_max = ddc::transform_reduce( + x_domain, + 0., + ddc::reducer::max(), + [](ddc::DiscreteElement ix) { + return 1. + / (ddc::distance_at_left(ix) + * ddc::distance_at_right(ix)); + }); + + double const invdy2_max = ddc::transform_reduce( + y_domain, + 0., + ddc::reducer::max(), + [](ddc::DiscreteElement iy) { + return 1. + / (ddc::distance_at_left(iy) + * ddc::distance_at_right(iy)); + }); + + ddc::Coordinate const max_dt { + .5 / (kx * invdx2_max + ky * invdy2_max)}; + + //! [CFL-condition] + + //! [time-domain] + ddc::DiscreteVector const nb_time_steps( + std::ceil((end_time - start_time) / max_dt) + .2); + + ddc::DiscreteDomain const time_domain + = ddc::init_discrete_space(DDimT::init( + ddc::Coordinate(start_time), + ddc::Coordinate(end_time), + nb_time_steps + 1)); + + //! [time-domain] + + //! [data allocation] + ddc::Chunk ghosted_last_temp( + "ghosted_last_temp", + ddc::DiscreteDomain< + DDimX, + DDimY>(ghosted_x_domain, ghosted_y_domain), + ddc::DeviceAllocator()); + + ddc::Chunk ghosted_next_temp( + "ghosted_next_temp", + ddc::DiscreteDomain< + DDimX, + DDimY>(ghosted_x_domain, ghosted_y_domain), + ddc::DeviceAllocator()); + //! [data allocation] + + //! [initial-chunkspan] + ddc::ChunkSpan const ghosted_initial_temp + = ghosted_last_temp.span_view(); + //! [initial-chunkspan] + + //! [fill-initial-chunkspan] + ddc::parallel_for_each( + ddc::DiscreteDomain(x_domain, y_domain), + KOKKOS_LAMBDA(ddc::DiscreteElement const ixy) { + double const x = ddc::coordinate( + ddc::DiscreteElement(ixy)); + double const y = ddc::coordinate( + ddc::DiscreteElement(ixy)); + ghosted_initial_temp(ixy) + = 9.999 * ((x * x + y * y) < 0.25); + }); + //! [fill-initial-chunkspan] + + //! [host-chunk] + ddc::Chunk ghosted_temp( + "ghost_temp", + ddc::DiscreteDomain< + DDimX, + DDimY>(ghosted_x_domain, ghosted_y_domain), + ddc::HostAllocator()); + //! [host-chunk] + + //! [initial-deepcopy] + ddc::parallel_deepcopy(ghosted_temp, ghosted_last_temp); + //! [initial-deepcopy] + + //! [initial-display] + display(ddc::coordinate(time_domain.front()), + ghosted_temp[x_domain][y_domain]); + //! [initial-display] + + //! [last-output-iter] + ddc::DiscreteElement last_output_iter = time_domain.front(); + //! [last-output-iter] + + //! [time iteration] + for (ddc::DiscreteElement const iter : + time_domain.remove_first(ddc::DiscreteVector(1))) { + //! [time iteration] + + //! [boundary conditions] + ddc::parallel_deepcopy( + ghosted_last_temp[ddc::DiscreteDomain< + DDimX, + DDimY>(x_pre_ghost, y_domain)], + ghosted_last_temp[ddc::DiscreteDomain< + DDimX, + DDimY>(y_domain, x_domain_end)]); + ddc::parallel_deepcopy( + ghosted_last_temp[ddc::DiscreteDomain< + DDimX, + DDimY>(y_domain, x_post_ghost)], + ghosted_last_temp[ddc::DiscreteDomain< + DDimX, + DDimY>(y_domain, x_domain_begin)]); + ddc::parallel_deepcopy( + ghosted_last_temp[ddc::DiscreteDomain< + DDimX, + DDimY>(x_domain, y_pre_ghost)], + ghosted_last_temp[ddc::DiscreteDomain< + DDimX, + DDimY>(x_domain, y_domain_end)]); + ddc::parallel_deepcopy( + ghosted_last_temp[ddc::DiscreteDomain< + DDimX, + DDimY>(x_domain, y_post_ghost)], + ghosted_last_temp[ddc::DiscreteDomain< + DDimX, + DDimY>(x_domain, y_domain_begin)]); + //! [boundary conditions] + + //! [manipulated views] + ddc::ChunkSpan const next_temp( + ghosted_next_temp[ddc::DiscreteDomain< + DDimX, + DDimY>(x_domain, y_domain)]); + ddc::ChunkSpan const last_temp(ghosted_last_temp.span_cview()); + //! [manipulated views] + + //! [numerical scheme] + ddc::parallel_for_each( + next_temp.domain(), + KOKKOS_LAMBDA( + ddc::DiscreteElement const ixy) { + ddc::DiscreteElement const ix + = ddc::select(ixy); + ddc::DiscreteElement const iy + = ddc::select(ixy); + double const dx_l = ddc::distance_at_left(ix); + double const dx_r = ddc::distance_at_right(ix); + double const dx_m = 0.5 * (dx_l + dx_r); + double const dy_l = ddc::distance_at_left(iy); + double const dy_r = ddc::distance_at_right(iy); + double const dy_m = 0.5 * (dy_l + dy_r); + next_temp(ix, iy) = last_temp(ix, iy); + next_temp(ix, iy) + += kx * ddc::step() + * (dx_l * last_temp(ix + 1, iy) + - 2.0 * dx_m * last_temp(ix, iy) + + dx_r * last_temp(ix - 1, iy)) + / (dx_l * dx_m * dx_r); + next_temp(ix, iy) + += ky * ddc::step() + * (dy_l * last_temp(ix, iy + 1) + - 2.0 * dy_m * last_temp(ix, iy) + + dy_r * last_temp(ix, iy - 1)) + / (dy_l * dy_m * dy_r); + }); + //! [numerical scheme] + + //! [output] + if (iter - last_output_iter >= t_output_period) { + last_output_iter = iter; + ddc::parallel_deepcopy(ghosted_temp, ghosted_last_temp); + display(ddc::coordinate(iter), + ghosted_temp[ddc::DiscreteDomain< + DDimX, + DDimY>(x_domain, y_domain)]); + } + //! [output] + + //! [swap] + std::swap(ghosted_last_temp, ghosted_next_temp); + //! [swap] + } + + + //! [final output] + if (last_output_iter < time_domain.back()) { + ddc::parallel_deepcopy(ghosted_temp, ghosted_last_temp); + display(ddc::coordinate(time_domain.back()), + ghosted_temp[ddc::DiscreteDomain< + DDimX, + DDimY>(x_domain, y_domain)]); + } + //! [final output] + +#ifdef DDC_BUILD_PDI_WRAPPER + PDI_finalize(); + PC_tree_destroy(&pdi_conf); +#endif +} diff --git a/examples/uniform_heat_equation.cpp b/examples/uniform_heat_equation.cpp index 91ea4623e..53c20d649 100644 --- a/examples/uniform_heat_equation.cpp +++ b/examples/uniform_heat_equation.cpp @@ -36,13 +36,17 @@ struct DDimT : ddc::UniformPointSampling }; //! [time-space] +//! [display] + /** A function to pretty print the temperature * @tparam ChunkType The type of chunk span. This way the template parameters are avoided, * should be deduced by the compiler. * @param time The time at which the output is made. * @param temp The temperature at this time-step. */ + //! [display] + template void display(double time, ChunkType temp) { @@ -107,6 +111,7 @@ int main(int argc, char** argv) //! [main-start-t-parameters] std::ptrdiff_t const t_output_period = 10; //! [parameters] + //! [main-start] //! [X-parameters] @@ -149,12 +154,14 @@ int main(int argc, char** argv) //! [Y-domains] //! [CFL-condition] + double const dx = ddc::step(); double const dy = ddc::step(); double const invdx2 = 1. / (dx * dx); double const invdy2 = 1. / (dy * dy); ddc::Coordinate const dt(.5 / (kx * invdx2 + ky * invdy2)); + //! [CFL-condition] //! [time-domain] @@ -166,6 +173,7 @@ int main(int argc, char** argv) ddc::Coordinate(start_time), ddc::Coordinate(end_time), nb_time_steps + 1)); + //! [time-domain] //! [data allocation] diff --git a/include/ddc/non_uniform_point_sampling.hpp b/include/ddc/non_uniform_point_sampling.hpp index 7bb957b4e..7a97d6708 100644 --- a/include/ddc/non_uniform_point_sampling.hpp +++ b/include/ddc/non_uniform_point_sampling.hpp @@ -113,6 +113,12 @@ class NonUniformPointSampling : detail::NonUniformPointSamplingBase return m_points.size(); } + /// @brief Lower bound index of the mesh + KOKKOS_FUNCTION discrete_element_type front() const noexcept + { + return discrete_element_type {0}; + } + /// @brief Convert a mesh index into a position in `CDim` KOKKOS_FUNCTION continuous_element_type coordinate(discrete_element_type const& icoord) const noexcept @@ -120,6 +126,78 @@ class NonUniformPointSampling : detail::NonUniformPointSamplingBase return m_points(icoord.uid()); } }; + + /** Construct an Impl and associated discrete_domain_type from an iterator + * containing the points coordinates along the `DDim` dimension. + * + * @param non_uniform_points a vector containing the coordinates of the points of the domain. + */ + template + static std::tuple, DiscreteDomain> + init(InputRange const non_uniform_points) + { + auto a = non_uniform_points.begin(); + auto b = non_uniform_points.end(); + auto n = std::distance(non_uniform_points.begin(), non_uniform_points.end()); + assert(a < b); + assert(n > 0); + typename DDim::template Impl disc(non_uniform_points); + DiscreteDomain domain {disc.front(), DiscreteVector {n}}; + return std::make_tuple(std::move(disc), std::move(domain)); + } + + /** Construct 4 non-uniform `DiscreteDomain` and an Impl from 3 iterators containing the points coordinates along the `DDim` dimension. + * + * @param domain_r an iterator containing the coordinates of the points of the main domain along the DDim position + * @param pre_ghost_r an iterator containing the positions of the ghost points before the main domain the DDim position + * @param post_ghost_r an iterator containing the positions of the ghost points after the main domain the DDim position + */ + template + static std::tuple< + typename DDim::template Impl, + DiscreteDomain, + DiscreteDomain, + DiscreteDomain, + DiscreteDomain> + init_ghosted( + InputRange const& domain_r, + InputRange const& pre_ghost_r, + InputRange const& post_ghost_r) + { + using discrete_domain_type = DiscreteDomain; + auto n = DiscreteVector {std::distance(domain_r.begin(), domain_r.end())}; + + assert(domain_r.begin() < domain_r.end()); + assert(n > 1); + + auto n_ghosts_before + = DiscreteVector {std::distance(pre_ghost_r.begin(), pre_ghost_r.end())}; + auto n_ghosts_after + = DiscreteVector {std::distance(post_ghost_r.begin(), post_ghost_r.end())}; + + std::vector full_domain; + + std::copy(pre_ghost_r.begin(), pre_ghost_r.end(), std::back_inserter(full_domain)); + std::copy(domain_r.begin(), domain_r.end(), std::back_inserter(full_domain)); + std::copy(post_ghost_r.begin(), post_ghost_r.end(), std::back_inserter(full_domain)); + + typename DDim::template Impl disc(full_domain); + + discrete_domain_type ghosted_domain + = discrete_domain_type(disc.front(), n + n_ghosts_before + n_ghosts_after); + discrete_domain_type pre_ghost + = discrete_domain_type(ghosted_domain.front(), n_ghosts_before); + discrete_domain_type main_domain + = discrete_domain_type(ghosted_domain.front() + n_ghosts_before, n); + discrete_domain_type post_ghost + = discrete_domain_type(main_domain.back() + 1, n_ghosts_after); + return std::make_tuple( + std::move(disc), + std::move(main_domain), + std::move(ghosted_domain), + std::move(pre_ghost), + std::move(post_ghost)); + } }; template From 15a155e11a63c076d46c44e481126b6b1dfd6c40 Mon Sep 17 00:00:00 2001 From: Thomas Padioleau Date: Mon, 24 Jun 2024 14:18:44 +0200 Subject: [PATCH 06/23] Update pages.yml --- .github/workflows/pages.yml | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/.github/workflows/pages.yml b/.github/workflows/pages.yml index e88213ba6..7d11084b1 100644 --- a/.github/workflows/pages.yml +++ b/.github/workflows/pages.yml @@ -161,7 +161,9 @@ jobs: mkdir docs_out chmod a+rwx docs_out docker run -v ${PWD}:/src ghcr.io/cexa-project/ddc/doxygen:${GITHUB_SHA:0:7} bash /src/run.sh - cat docs_out/doxygen.log + if [ -f docs_out/doxygen.log ]; then + cat docs_out/doxygen.log + fi - name: Publish site if: ${{ github.event_name == 'push' && github.ref_name == 'main' && needs.id_repo.outputs.in_base_repo == 'true' }} run: | From 5fd96f85f52690b16482296767ba71da074b8233 Mon Sep 17 00:00:00 2001 From: Baptiste Legouix Date: Tue, 25 Jun 2024 12:57:24 +0200 Subject: [PATCH 07/23] Splines minor classes doc (#452) --- .../ddc/kernels/splines/greville_interpolation_points.hpp | 6 ++++++ .../ddc/kernels/splines/knots_as_interpolation_points.hpp | 4 ++++ 2 files changed, 10 insertions(+) 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 From 6b4b5b96fc55c47cb8625bbccf799ffe4a16089b Mon Sep 17 00:00:00 2001 From: Baptiste Legouix Date: Tue, 25 Jun 2024 18:17:58 +0200 Subject: [PATCH 08/23] SplineEvaluator doc (#430) * init * missing sed in benchmark * fix * Update include/ddc/kernels/splines/spline_builder.hpp * clang-format * fix kokkos version change * wip * improve consistency * wip * wip * spline_builder * changes from Emiliy's reviw * format * minor * wip * wip * wip * bsplines * init * reinit * remove bsplines (not in the scope of this MR anymore) * format * autoreview * wip * wip * wip * fix for allowing doxygen to keep track of integrals() * format * Revert "fix for allowing doxygen to keep track of integrals()" This reverts commit 70172346a26721f023620e04ad15a60d2296985e. * autoreview * B-splines * Update include/ddc/kernels/splines/bsplines_non_uniform.hpp Co-authored-by: EmilyBourne * emily's review * wip * Update include/ddc/kernels/splines/bsplines_non_uniform.hpp Co-authored-by: Thomas Padioleau * wip * reviews * typo * fix * minor * more details on uniformity * non-uniform constructors * minor * non_uniform constructors again * wip * wip on CI * wip on CI * CI, integrals doxygen tracking still wrong * fix * Update include/ddc/kernels/splines/bsplines_non_uniform.hpp Co-authored-by: EmilyBourne * wip * emily's review * minor * minor * the spline coefficients * CI * emily's review * forgot a files save * fix doxygen * clang-format * wip on CI * null_extrap * CI * wip * Emily's review * wip * Thomas' review * autoreview * init * remove null_extrapolation * Emily's minireview * shorten doxygen comments * wip * Emily's review * ident * wip * wip * wip * wip * wip * CI * wip * Emily's review * minor * ident * ident * minor * transform -> approximate * wip * Emily's review * init * caps * minor * Emily's review * Update include/ddc/kernels/splines/spline_builder.hpp Co-authored-by: EmilyBourne * Emily's review * minor * Emily's review * minor * interest coordinate * wip * wip * wip * SplineEvaluator2D * autoreview * ident * hyperparameter -> parameter * init * caps * interest coordinate * Apply suggestions from code review * Apply suggestions from code review Co-authored-by: Thomas Padioleau --------- Co-authored-by: EmilyBourne Co-authored-by: Thomas Padioleau --- .../splines/null_extrapolation_rule.hpp | 8 + .../ddc/kernels/splines/spline_evaluator.hpp | 175 +++++++- .../kernels/splines/spline_evaluator_2d.hpp | 397 +++++++++++++----- 3 files changed, 472 insertions(+), 108 deletions(-) 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_evaluator.hpp b/include/ddc/kernels/splines/spline_evaluator.hpp index d85d64b1b..444116498 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 @@ -170,6 +287,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 +307,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 @@ -208,6 +354,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, diff --git a/include/ddc/kernels/splines/spline_evaluator_2d.hpp b/include/ddc/kernels/splines/spline_evaluator_2d.hpp index fa77349a2..57701ead2 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 @@ -304,14 +395,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 +415,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 +435,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 +454,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 +489,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 +524,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( @@ -443,14 +576,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( @@ -485,14 +626,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( @@ -526,6 +675,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 +722,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 +775,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( From d6e6b58d3e33e7ff687b026c2e6b45bf0b495bb3 Mon Sep 17 00:00:00 2001 From: Baptiste Legouix Date: Wed, 26 Jun 2024 11:12:49 +0200 Subject: [PATCH 09/23] Splines linear problem periodic band (#477) --- .../splines/splines_linear_problem_maker.hpp | 35 +++++++++++++++++++ tests/splines/matrix.cpp | 29 +++++++++++++++ 2 files changed, 64 insertions(+) diff --git a/include/ddc/kernels/splines/splines_linear_problem_maker.hpp b/include/ddc/kernels/splines/splines_linear_problem_maker.hpp index 3823ac89c..2826f5692 100644 --- a/include/ddc/kernels/splines/splines_linear_problem_maker.hpp +++ b/include/ddc/kernels/splines/splines_linear_problem_maker.hpp @@ -93,6 +93,41 @@ class SplinesLinearProblemMaker SplinesLinearProblem2x2Blocks>(n, std::move(top_left_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); + } + /** * @brief Construct a sparse matrix * diff --git a/tests/splines/matrix.cpp b/tests/splines/matrix.cpp index a1e8b1fec..e99d5d78e 100644 --- a/tests/splines/matrix.cpp +++ b/tests/splines/matrix.cpp @@ -339,6 +339,35 @@ TEST_P(MatrixSizesFixture, 2x2Blocks) 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(); From 6071a80a04dedc8fa9a4c61beae27093a4968dc3 Mon Sep 17 00:00:00 2001 From: thierryantoun <152993060+thierryantoun@users.noreply.github.com> Date: Wed, 26 Jun 2024 16:04:19 +0200 Subject: [PATCH 10/23] type correction in the non uniform heat equations (#505) --- examples/non_uniform_heat_equation.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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()); From 46c6f5e35ae45888c9b5b638655e8b04be88c1a5 Mon Sep 17 00:00:00 2001 From: Baptiste Legouix Date: Wed, 26 Jun 2024 17:37:32 +0200 Subject: [PATCH 11/23] Apply suggestions from code review Co-authored-by: Thomas Padioleau --- .../splines/splines_linear_problem_3x3_blocks.hpp | 14 +++----------- 1 file changed, 3 insertions(+), 11 deletions(-) diff --git a/include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp b/include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp index 8c4d0f202..977c1f81f 100644 --- a/include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp +++ b/include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp @@ -95,12 +95,12 @@ class SplinesLinearProblem3x3Blocks : public SplinesLinearProblem2x2Blocks | b_top | -- Considered as a + * | 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 b) const + 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 @@ -110,10 +110,6 @@ class SplinesLinearProblem3x3Blocks : public SplinesLinearProblem2x2Blocks {m_top_size, m_top_size + nq}, Kokkos::ALL); - MultiRHS const b_bottom = Kokkos:: - subview(b, - std::pair {m_top_size + nq, size()}, - Kokkos::ALL); MultiRHS const b_center_dst = Kokkos::subview(b, std::pair {0, nq}, Kokkos::ALL); @@ -151,10 +147,6 @@ class SplinesLinearProblem3x3Blocks : public SplinesLinearProblem2x2Blocks {m_top_size, m_top_size + nq}, Kokkos::ALL); - MultiRHS const b_bottom = Kokkos:: - subview(b, - std::pair {m_top_size + nq, size()}, - Kokkos::ALL); MultiRHS const buffer = Kokkos::create_mirror(ExecSpace(), b_center); @@ -172,7 +164,7 @@ class SplinesLinearProblem3x3Blocks : public SplinesLinearProblem2x2Blocks Date: Wed, 26 Jun 2024 17:37:44 +0200 Subject: [PATCH 12/23] Update include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp Co-authored-by: Thomas Padioleau --- .../ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp b/include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp index 977c1f81f..de8bb0e52 100644 --- a/include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp +++ b/include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp @@ -132,7 +132,7 @@ class SplinesLinearProblem3x3Blocks : public SplinesLinearProblem2x2Blockssize(); // size of the center block From 3c2da245e1b745d5dd42f326e09604f32e042ab4 Mon Sep 17 00:00:00 2001 From: Baptiste Legouix Date: Wed, 26 Jun 2024 17:38:07 +0200 Subject: [PATCH 13/23] Apply suggestions from code review Co-authored-by: Thomas Padioleau --- .../ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp b/include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp index de8bb0e52..e7a7d75cc 100644 --- a/include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp +++ b/include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp @@ -159,7 +159,7 @@ class SplinesLinearProblem3x3Blocks : public SplinesLinearProblem2x2Blocks Date: Wed, 26 Jun 2024 17:38:40 +0200 Subject: [PATCH 14/23] Update include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp Co-authored-by: Thomas Padioleau --- .../splines/splines_linear_problem_3x3_blocks.hpp | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp b/include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp index e7a7d75cc..c2bcdb621 100644 --- a/include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp +++ b/include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp @@ -65,15 +65,17 @@ class SplinesLinearProblem3x3Blocks : public SplinesLinearProblem2x2Blockssize(); // size of the center block - if (i < m_top_size) + if (i < m_top_size) { i += nq; - else if (i < m_top_size + nq) + } else if (i < m_top_size + nq) { i -= m_top_size; + } - if (j < m_top_size) + if (j < m_top_size) { j += nq; - else if (j < m_top_size + nq) + } else if (j < m_top_size + nq) { j -= m_top_size; + } } public: From c56c130009a6e0a1839f97d0acdce2eb1ea960c2 Mon Sep 17 00:00:00 2001 From: blegouix Date: Wed, 26 Jun 2024 17:57:30 +0200 Subject: [PATCH 15/23] asserts to prevent overlapping allocations --- .../kernels/splines/splines_linear_problem_3x3_blocks.hpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp b/include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp index c2bcdb621..1dc93e99b 100644 --- a/include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp +++ b/include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp @@ -106,6 +106,9 @@ class SplinesLinearProblem3x3Blocks : public SplinesLinearProblem2x2Blockssize(); // 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:: @@ -138,6 +141,9 @@ class SplinesLinearProblem3x3Blocks : public SplinesLinearProblem2x2Blockssize(); // 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:: From 93a2b2d3644a09f5668a34a6cb38c6be9549f62f Mon Sep 17 00:00:00 2001 From: blegouix Date: Wed, 26 Jun 2024 18:07:26 +0200 Subject: [PATCH 16/23] reset 2x2 file --- .../splines_linear_problem_2x2_blocks.hpp | 32 ++----------------- 1 file changed, 2 insertions(+), 30 deletions(-) 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 3769c5f07..c459166fb 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 it. + * @param top_left_block A pointer toward the top-left SplinesLinearProblem. `setup_solver` must not have been called on `q`. */ explicit SplinesLinearProblem2x2Blocks( std::size_t const mat_size, @@ -74,34 +74,6 @@ class SplinesLinearProblem2x2Blocks : public SplinesLinearProblem Kokkos::deep_copy(m_bottom_left_block.h_view, 0.); } -protected: - /** - * @brief SplinesLinearProblem2x2Blocks constructor. - * - * @param mat_size The size of one of the dimensions of the square matrix. - * @param q A pointer toward the top-left SplinesLinearProblem. - */ - explicit SplinesLinearProblem2x2Blocks( - std::size_t const mat_size, - std::unique_ptr> top_left_block, - std::size_t const lambda_size1, - std::size_t const lambda_size2) - : SplinesLinearProblem(mat_size) - , m_top_left_block(std::move(top_left_block)) - , m_top_right_block( - "top_right_block", - m_top_left_block->size(), - mat_size - m_top_left_block->size()) - , m_bottom_left_block("bottom_left_block", lambda_size1, lambda_size2) - , m_bottom_right_block( - new SplinesLinearProblemDense(mat_size - m_top_left_block->size())) - { - assert(m_top_left_block->size() <= mat_size); - - Kokkos::deep_copy(m_top_right_block.h_view, 0.); - Kokkos::deep_copy(m_bottom_left_block.h_view, 0.); - } - double get_element(std::size_t const i, std::size_t const j) const override { assert(i < size()); @@ -201,7 +173,7 @@ class SplinesLinearProblem2x2Blocks : public SplinesLinearProblem * @param LinOp * @param transpose */ - virtual void gemv_minus1_1( + void gemv_minus1_1( MultiRHS const x, MultiRHS const y, Kokkos::View const From 20ea4696f8b4795b7aa2f1deca0fd1a0c085f6d4 Mon Sep 17 00:00:00 2001 From: Baptiste Legouix Date: Wed, 26 Jun 2024 18:58:06 +0200 Subject: [PATCH 17/23] Update include/ddc/kernels/splines/splines_linear_problem_2x2_blocks.hpp --- .../ddc/kernels/splines/splines_linear_problem_2x2_blocks.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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..7b7461e99 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, From e3abfe8e9264d73064b845fbfd4a755fa2a839ff Mon Sep 17 00:00:00 2001 From: Baptiste Legouix Date: Wed, 26 Jun 2024 18:58:29 +0200 Subject: [PATCH 18/23] Update include/ddc/kernels/splines/splines_linear_problem_2x2_blocks.hpp --- .../ddc/kernels/splines/splines_linear_problem_2x2_blocks.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 7b7461e99..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 `it`. + * @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, From 86484b7599b8a970a0b523598fc8982c68188c8e Mon Sep 17 00:00:00 2001 From: Baptiste Legouix Date: Wed, 26 Jun 2024 18:58:58 +0200 Subject: [PATCH 19/23] Update tests/splines/CMakeLists.txt --- tests/splines/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/splines/CMakeLists.txt b/tests/splines/CMakeLists.txt index f707a22f2..54b717828 100644 --- a/tests/splines/CMakeLists.txt +++ b/tests/splines/CMakeLists.txt @@ -113,7 +113,7 @@ foreach(BC "BC_PERIODIC" "BC_GREVILLE" "BC_HERMITE") GTest::gtest DDC::DDC ) - target_compile_definitions("${test_name}" PUBLIC -DDEGREE_X=${DEGREE_X} -D${BSPLINES_TYPE} -D${BC} -DSOLVER_LAPACK) + target_compile_definitions("${test_name}" PUBLIC -DDEGREE_X=${DEGREE_X} -D${BSPLINES_TYPE} -D${BC} -DSOLVER_LAPACK) # add_test("${test_name}" "${test_name}") gtest_discover_tests("${test_name}" DISCOVERY_MODE PRE_TEST) endforeach() From 8f8dc8dd6aaf4f8ca896a08417759171cc6c1d4d Mon Sep 17 00:00:00 2001 From: Baptiste Legouix Date: Wed, 26 Jun 2024 22:09:34 +0200 Subject: [PATCH 20/23] Splines linear problem 3x3 blocks (#478) * wip * wip * runs but wrong result * fix * rely on Kokkos::View to store nb and k * renaming * misc * release candidate * autoreview * autoreview * wip * wip * thomas' review * improve gemv * minor fix * minor fix * renaming * fix * init * autoreview + missing files * autoreview * autoreview * Apply suggestions from code review Co-authored-by: Thomas Padioleau * Update include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp Co-authored-by: Thomas Padioleau * Apply suggestions from code review Co-authored-by: Thomas Padioleau * Update include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp Co-authored-by: Thomas Padioleau * asserts to prevent overlapping allocations --------- Co-authored-by: Thomas Padioleau --- include/ddc/kernels/splines.hpp | 1 + .../splines_linear_problem_2x2_blocks.hpp | 2 +- .../splines_linear_problem_3x3_blocks.hpp | 185 ++++++++++++++++++ .../splines/splines_linear_problem_maker.hpp | 23 ++- tests/splines/matrix.cpp | 48 ++++- 5 files changed, 249 insertions(+), 10 deletions(-) create mode 100644 include/ddc/kernels/splines/splines_linear_problem_3x3_blocks.hpp 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/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 2826f5692..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,18 @@ 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)); } /** diff --git a/tests/splines/matrix.cpp b/tests/splines/matrix.cpp index e99d5d78e..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,27 @@ 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(); From be38df7f072cdf4b9beaedddc3bfc8a41a71c077 Mon Sep 17 00:00:00 2001 From: Thomas Padioleau Date: Thu, 27 Jun 2024 07:01:06 +0200 Subject: [PATCH 21/23] Add readability-avoid-const-params-in-decls (#507) --- .clang-tidy | 2 +- .../ddc/kernels/splines/spline_builder.hpp | 16 ++--- .../ddc/kernels/splines/spline_builder_2d.hpp | 64 +++++++------------ 3 files changed, 31 insertions(+), 51 deletions(-) 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/include/ddc/kernels/splines/spline_builder.hpp b/include/ddc/kernels/splines/spline_builder.hpp index 9f0332e9c..502d099e7 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: 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; }; From 9e6ea8f49caf7b981ce73da5faaaea9f16264116 Mon Sep 17 00:00:00 2001 From: Baptiste Legouix Date: Thu, 27 Jun 2024 14:05:09 +0200 Subject: [PATCH 22/23] Labelize ddc kernels (#509) --- include/ddc/kernels/fft.hpp | 1 + include/ddc/kernels/splines/spline_builder.hpp | 12 +++++++++--- include/ddc/kernels/splines/spline_evaluator.hpp | 4 ++++ include/ddc/kernels/splines/spline_evaluator_2d.hpp | 6 ++++++ 4 files changed, 20 insertions(+), 3 deletions(-) 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/spline_builder.hpp b/include/ddc/kernels/splines/spline_builder.hpp index 502d099e7..697deb5d3 100644 --- a/include/ddc/kernels/splines/spline_builder.hpp +++ b/include/ddc/kernels/splines/spline_builder.hpp @@ -766,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) { @@ -783,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) { @@ -803,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) { @@ -816,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) { @@ -837,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) { @@ -848,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_evaluator.hpp b/include/ddc/kernels/splines/spline_evaluator.hpp index 444116498..97458707a 100644 --- a/include/ddc/kernels/splines/spline_evaluator.hpp +++ b/include/ddc/kernels/splines/spline_evaluator.hpp @@ -275,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) { @@ -341,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) { @@ -379,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 57701ead2..53c13b74c 100644 --- a/include/ddc/kernels/splines/spline_evaluator_2d.hpp +++ b/include/ddc/kernels/splines/spline_evaluator_2d.hpp @@ -380,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) { @@ -559,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) { @@ -609,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) { @@ -659,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) { @@ -804,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); @@ -811,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) { From 72fc3d6af05f7d33fad2785808d3c1535d810328 Mon Sep 17 00:00:00 2001 From: blegouix Date: Fri, 28 Jun 2024 11:07:56 +0200 Subject: [PATCH 23/23] restore ginkgo tests --- tests/splines/CMakeLists.txt | 4 ++-- tests/splines/batched_spline_builder.cpp | 10 ++++++++++ 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/tests/splines/CMakeLists.txt b/tests/splines/CMakeLists.txt index 54b717828..de6476f24 100644 --- a/tests/splines/CMakeLists.txt +++ b/tests/splines/CMakeLists.txt @@ -121,9 +121,9 @@ foreach(BC "BC_PERIODIC" "BC_GREVILLE" "BC_HERMITE") endforeach() # GINKGO -foreach(BC "BC_PERIODIC") +foreach(BC "BC_PERIODIC" "BC_GREVILLE" "BC_HERMITE") foreach(DEGREE_X RANGE "${SPLINES_TEST_DEGREE_MIN}" "${SPLINES_TEST_DEGREE_MAX}") - foreach(BSPLINES_TYPE "BSPLINES_TYPE_UNIFORM") + foreach(BSPLINES_TYPE "BSPLINES_TYPE_UNIFORM" "BSPLINES_TYPE_NON_UNIFORM") set(test_name "splines_tests_BATCHED_DEGREE_X_${DEGREE_X}_${BSPLINES_TYPE}_${BC}_GINKGO") add_executable("${test_name}" ../main.cpp batched_spline_builder.cpp) target_compile_features("${test_name}" PUBLIC cxx_std_17) diff --git a/tests/splines/batched_spline_builder.cpp b/tests/splines/batched_spline_builder.cpp index c2d57805e..41d968a60 100644 --- a/tests/splines/batched_spline_builder.cpp +++ b/tests/splines/batched_spline_builder.cpp @@ -434,6 +434,16 @@ static void BatchedSplineTest() #define SUFFIX(name) name##Lapack##Hermite##NonUniform #elif defined(BC_PERIODIC) && defined(BSPLINES_TYPE_UNIFORM) && defined(SOLVER_GINKGO) #define SUFFIX(name) name##Ginkgo##Periodic##Uniform +#elif defined(BC_PERIODIC) && defined(BSPLINES_TYPE_NON_UNIFORM) && defined(SOLVER_GINKGO) +#define SUFFIX(name) name##Ginkgo##Periodic##NonUniform +#elif defined(BC_GREVILLE) && defined(BSPLINES_TYPE_UNIFORM) && defined(SOLVER_GINKGO) +#define SUFFIX(name) name##Ginkgo##Greville##Uniform +#elif defined(BC_GREVILLE) && defined(BSPLINES_TYPE_NON_UNIFORM) && defined(SOLVER_GINKGO) +#define SUFFIX(name) name##Ginkgo##Greville##NonUniform +#elif defined(BC_HERMITE) && defined(BSPLINES_TYPE_UNIFORM) && defined(SOLVER_GINKGO) +#define SUFFIX(name) name##Ginkgo##Hermite##Uniform +#elif defined(BC_HERMITE) && defined(BSPLINES_TYPE_NON_UNIFORM) && defined(SOLVER_GINKGO) +#define SUFFIX(name) name##Ginkgo##Hermite##NonUniform #endif TEST(SUFFIX(BatchedSplineHost), 1DX)