wip

CExA-project · Jul 23, 2024 · 03f45a3 · 03f45a3
1 parent eb8ed5f
commit 03f45a3
Show file tree

Hide file tree

Showing 2 changed files with 39 additions and 17 deletions.
diff --git a/include/ddc/kernels/splines/spline_builder.hpp b/include/ddc/kernels/splines/spline_builder.hpp
@@ -401,8 +401,7 @@ class SplineBuilder
             double,
             ddc::DiscreteDomain<interpolation_discrete_dimension_type>,
             ddc::KokkosAllocator<double, memory_space>>
-    quadrature_coefficients(
-            ddc::DiscreteDomain<interpolation_discrete_dimension_type> const& domain) const;
+    quadrature_coefficients() const;
 
 private:
     void compute_block_sizes_uniform(int& lower_block_size, int& upper_block_size) const;
@@ -922,17 +921,19 @@ SplineBuilder<
         BcLower,
         BcUpper,
         Solver,
-        IDimX...>::quadrature_coefficients(ddc::DiscreteDomain<InterpolationDDim> const& domain)
-        const
+        IDimX...>::quadrature_coefficients() const
 {
+    // Compute integrals of bsplines
     ddc::Chunk<double, ddc::DiscreteDomain<bsplines_type>> integral_bsplines(spline_domain());
     ddc::discrete_space<bsplines_type>().integrals(integral_bsplines.span_view());
 
-    // Solve matrix equation
+    // Remove last cell in the periodic case
     ddc::ChunkSpan integral_bsplines_without_periodic_point
             = integral_bsplines.span_view()[ddc::DiscreteDomain<bsplines_type>(
                     ddc::DiscreteElement<bsplines_type>(0),
                     ddc::DiscreteVector<bsplines_type>(matrix->size()))];
+
+    // Allocate
     Kokkos::View<double**, Kokkos::LayoutRight, Kokkos::DefaultHostExecutionSpace>
             integral_bsplines_mirror_with_additional_allocation(
                     "integral_bsplines_mirror_with_additional_allocation",
@@ -945,6 +946,8 @@ SplineBuilder<
                                     0,
                                     integral_bsplines_without_periodic_point.size()},
                             0);
+
+    // Solve matrix equation A^t*X=integral_bsplines
     Kokkos::deep_copy(
             integral_bsplines_mirror,
             integral_bsplines_without_periodic_point.allocation_kokkos_view());
@@ -957,18 +960,17 @@ SplineBuilder<
             double,
             ddc::DiscreteDomain<InterpolationDDim>,
             ddc::KokkosAllocator<double, MemorySpace>>
-            coefficients(domain);
-
-    // Coefficients of quadrature in integral_bsplines (values which would always be multiplied
-    // by f'(x)=0 are removed
-    ddc::DiscreteDomain<bsplines_type> slice =
-            spline_domain()
-                      .remove(ddc::DiscreteVector<bsplines_type> {s_nbc_xmin},
-                              ddc::DiscreteVector<bsplines_type> {s_nbc_xmax});
-
+            coefficients(ddc::DiscreteDomain<interpolation_discrete_dimension_type>(
+                    interpolation_domain().front(),
+                    ddc::DiscreteVector<interpolation_discrete_dimension_type>(
+                            integral_bsplines_without_periodic_point.size())));
     Kokkos::deep_copy(
             coefficients.allocation_kokkos_view(),
-            integral_bsplines_without_periodic_point[slice].allocation_kokkos_view());
+            integral_bsplines_without_periodic_point
+                    [ddc::DiscreteDomain<bsplines_type>(
+                             ddc::DiscreteElement<bsplines_type>(s_nbc_xmin),
+                             ddc::DiscreteVector<bsplines_type>(interpolation_domain().size()))]
+                            .allocation_kokkos_view());
 
     return coefficients;
 }

diff --git a/tests/splines/periodic_spline_builder.cpp b/tests/splines/periodic_spline_builder.cpp
@@ -129,6 +129,21 @@ TEST(PeriodicSplineBuilderTest, Identity)
     spline_evaluator
             .deriv(spline_eval_deriv.span_view(), coords_eval.span_cview(), coef.span_cview());
 
+    ddc::Chunk integral(spline_builder.batch_domain(), ddc::HostAllocator<double>());
+    spline_evaluator.integrate(integral.span_view(), coef.span_cview());
+
+    ddc::Chunk quadrature_coefficients = spline_builder.quadrature_coefficients();
+    ddc::Chunk quadrature_integral(spline_builder.batch_domain(), ddc::HostAllocator<double>());
+    quadrature_integral(ddc::DiscreteElement<>()) = ddc::parallel_transform_reduce(
+            Kokkos::DefaultHostExecutionSpace(),
+            quadrature_coefficients.domain(),
+            0.0,
+            ddc::reducer::sum<double>(),
+            [&](ddc::DiscreteElement<IDimX> const ix) {
+                return quadrature_coefficients(ix) * yvals(ix);
+            });
+
+
     // 8. Checking errors
     std::cout << "---------- TEST ----------\n";
     double max_norm_error = 0.;
@@ -144,10 +159,12 @@ TEST(PeriodicSplineBuilderTest, Identity)
         double const error_deriv = spline_eval_deriv(ix) - evaluator.deriv(x, 1);
         max_norm_error_diff = std::fmax(max_norm_error_diff, std::fabs(error_deriv));
     }
-    ddc::Chunk integral(spline_builder.batch_domain(), ddc::HostAllocator<double>());
-    spline_evaluator.integrate(integral.span_view(), coef.span_cview());
+
     double const max_norm_error_integ = std::fabs(
             integral(ddc::DiscreteElement<>()) - evaluator.deriv(xN, -1) + evaluator.deriv(x0, -1));
+    double const max_norm_error_quadrature_integ = std::fabs(
+            quadrature_integral(ddc::DiscreteElement<>()) - evaluator.deriv(xN, -1)
+            + evaluator.deriv(x0, -1));
 
     double const max_norm = evaluator.max_norm();
     double const max_norm_diff = evaluator.max_norm(1);
@@ -164,4 +181,7 @@ TEST(PeriodicSplineBuilderTest, Identity)
     EXPECT_LE(
             max_norm_error_integ,
             std::max(error_bounds.error_bound_on_int(h, s_degree_x), 1.0e-14 * max_norm_int));
+    EXPECT_LE(
+            max_norm_error_quadrature_integ,
+            std::max(error_bounds.error_bound_on_int(h, s_degree_x), 1.0e-14 * max_norm_int));
 }