Provide a new initialization interface

CExA-project · Dec 29, 2024 · d67104a · d67104a
1 parent 2f8726a
commit d67104a
Show file tree

Hide file tree

Showing 33 changed files with 413 additions and 215 deletions.
diff --git a/benchmarks/splines.cpp b/benchmarks/splines.cpp
@@ -111,26 +111,24 @@ void characteristics_advection_unitary(benchmark::State& state)
             std::ref(maxUsedMem));
 
     if constexpr (!IsNonUniform) {
-        ddc::init_discrete_space<BSplinesX<
+        ddc::create_uniform_bsplines<BSplinesX<
                 IsNonUniform,
                 s_degree_x>>(ddc::Coordinate<X>(0.), ddc::Coordinate<X>(1.), nx);
     } else {
         std::vector<ddc::Coordinate<X>> breaks(nx + 1);
         for (std::size_t i(0); i < nx + 1; ++i) {
             breaks[i] = ddc::Coordinate<X>(static_cast<double>(i) / nx);
         }
-        ddc::init_discrete_space<BSplinesX<IsNonUniform, s_degree_x>>(breaks);
+        ddc::create_non_uniform_bsplines<BSplinesX<IsNonUniform, s_degree_x>>(breaks);
     }
-    ddc::init_discrete_space<DDimX<IsNonUniform, s_degree_x>>(
+    ddc::init_discrete_space_from_impl<DDimX<IsNonUniform, s_degree_x>>(
             ddc::GrevilleInterpolationPoints<
                     BSplinesX<IsNonUniform, s_degree_x>,
                     ddc::BoundCond::PERIODIC,
                     ddc::BoundCond::PERIODIC>::
                     template get_sampling<DDimX<IsNonUniform, s_degree_x>>());
-    ddc::DiscreteDomain<DDimY> const y_domain = ddc::init_discrete_space<DDimY>(DDimY::init<DDimY>(
-            ddc::Coordinate<Y>(-1.),
-            ddc::Coordinate<Y>(1.),
-            ddc::DiscreteVector<DDimY>(ny)));
+    ddc::DiscreteDomain<DDimY> const y_domain = ddc::create_uniform_point_sampling<
+            DDimY>(ddc::Coordinate<Y>(-1.), ddc::Coordinate<Y>(1.), ddc::DiscreteVector<DDimY>(ny));
 
     auto const x_domain = ddc::GrevilleInterpolationPoints<
             BSplinesX<IsNonUniform, s_degree_x>,

diff --git a/docs/first_steps.md b/docs/first_steps.md
@@ -75,7 +75,7 @@ Their type is not specified because we use C++
 all of the same type: `DiscreteDomain<DDimX>` that represents a set of contiguous points in the
 discretization of `X`.
 
-\ref ddc::UniformPointSampling::init_ghosted "init_ghosted" takes as parameters the coordinate of the first and last discretized points, the
+\ref ddc::create_uniform_point_sampling_ghosted "init_ghosted" takes as parameters the coordinate of the first and last discretized points, the
 number of discretized points in the domain and the number of additional points on each side of the
 domain.
 The fours `DiscreteDomain`s returned are:

diff --git a/examples/characteristics_advection.cpp b/examples/characteristics_advection.cpp
@@ -133,17 +133,17 @@ int main(int argc, char** argv)
 
     //! [X-global-domain]
     // Initialization of the global domain in X
-    ddc::init_discrete_space<
+    ddc::create_uniform_bsplines<
             BSplinesX>(ddc::Coordinate<X>(x_start), ddc::Coordinate<X>(x_end), nb_x_points);
-    ddc::init_discrete_space<DDimX>(GrevillePoints::get_sampling<DDimX>());
+    ddc::init_discrete_space_from_impl<DDimX>(GrevillePoints::get_sampling<DDimX>());
 
     auto const x_domain = GrevillePoints::get_domain<DDimX>();
     //! [X-global-domain]
     // Initialization of the global domain in Y
-    auto const y_domain = ddc::init_discrete_space<DDimY>(DDimY::init<DDimY>(
+    ddc::DiscreteDomain<DDimY> const y_domain = ddc::create_uniform_point_sampling<DDimY>(
             ddc::Coordinate<Y>(y_start),
             ddc::Coordinate<Y>(y_end),
-            ddc::DiscreteVector<DDimY>(nb_y_points)));
+            ddc::DiscreteVector<DDimY>(nb_y_points));
 
     //! [time-domains]
 
@@ -153,11 +153,8 @@ int main(int argc, char** argv)
     // Initialization of the global domain in time:
     // - the number of discrete time-points is equal to the number of
     //   steps + 1
-    ddc::DiscreteDomain<DDimT> const time_domain
-            = ddc::init_discrete_space<DDimT>(DDimT::init<DDimT>(
-                    ddc::Coordinate<T>(start_time),
-                    ddc::Coordinate<T>(end_time),
-                    nb_time_steps + 1));
+    ddc::DiscreteDomain<DDimT> const time_domain = ddc::create_uniform_point_sampling<
+            DDimT>(ddc::Coordinate<T>(start_time), ddc::Coordinate<T>(end_time), nb_time_steps + 1);
     //! [time-domains]
 
     //! [data allocation]

diff --git a/examples/heat_equation.cpp b/examples/heat_equation.cpp
@@ -122,11 +122,11 @@ int main(int argc, char** argv)
     // Initialization of the global domain in X with gwx ghost points on
     // each side
     auto const [x_domain, ghosted_x_domain, x_pre_ghost, x_post_ghost]
-            = ddc::init_discrete_space<DDimX>(DDimX::init_ghosted<DDimX>(
+            = ddc::create_uniform_point_sampling_ghosted<DDimX>(
                     ddc::Coordinate<X>(x_start),
                     ddc::Coordinate<X>(x_end),
                     ddc::DiscreteVector<DDimX>(nb_x_points),
-                    gwx));
+                    gwx);
     //! [X-global-domain]
 
     //! [X-domains]
@@ -147,11 +147,11 @@ int main(int argc, char** argv)
     // Initialization of the global domain in Y with gwy ghost points on
     // each side
     auto const [y_domain, ghosted_y_domain, y_pre_ghost, y_post_ghost]
-            = ddc::init_discrete_space<DDimY>(DDimY::init_ghosted<DDimY>(
+            = ddc::create_uniform_point_sampling_ghosted<DDimY>(
                     ddc::Coordinate<Y>(y_start),
                     ddc::Coordinate<Y>(y_end),
                     ddc::DiscreteVector<DDimY>(nb_y_points),
-                    gwy));
+                    gwy);
 
     // our zone at the start of the domain that will be mirrored to the
     // ghost
@@ -190,11 +190,8 @@ int main(int argc, char** argv)
     // - the number of discrete time-points is equal to the number of
     //   steps + 1
     // `init` takes required information to initialize the attributes of the dimension.
-    ddc::DiscreteDomain<DDimT> const time_domain
-            = ddc::init_discrete_space<DDimT>(DDimT::init<DDimT>(
-                    ddc::Coordinate<T>(start_time),
-                    ddc::Coordinate<T>(end_time),
-                    nb_time_steps + 1));
+    ddc::DiscreteDomain<DDimT> const time_domain = ddc::create_uniform_point_sampling<
+            DDimT>(ddc::Coordinate<T>(start_time), ddc::Coordinate<T>(end_time), nb_time_steps + 1);
     //! [time-domains]
 
     //! [data allocation]

diff --git a/examples/heat_equation_spectral.cpp b/examples/heat_equation_spectral.cpp
@@ -99,18 +99,18 @@ int main(int argc, char** argv)
     std::ptrdiff_t const t_output_period = 10;
 
     // Initialization of the global domain in X including periodic point to have correct step
-    auto const x_domain_with_periodic_point = ddc::init_discrete_space<DDimX>(DDimX::init<DDimX>(
+    auto const x_domain_with_periodic_point = ddc::create_uniform_point_sampling<DDimX>(
             ddc::Coordinate<X>(x_start),
             ddc::Coordinate<X>(x_end),
-            ddc::DiscreteVector<DDimX>(nb_x_points)));
+            ddc::DiscreteVector<DDimX>(nb_x_points));
     ddc::DiscreteDomain<DDimX> const x_domain
             = x_domain_with_periodic_point.remove_last(ddc::DiscreteVector<DDimX>(1));
 
     // Initialization of the global domain in Y including periodic point to have correct step
-    auto const y_domain_with_periodic_point = ddc::init_discrete_space<DDimY>(DDimY::init<DDimY>(
+    auto const y_domain_with_periodic_point = ddc::create_uniform_point_sampling<DDimY>(
             ddc::Coordinate<Y>(y_start),
             ddc::Coordinate<Y>(y_end),
-            ddc::DiscreteVector<DDimY>(nb_y_points)));
+            ddc::DiscreteVector<DDimY>(nb_y_points));
     ddc::DiscreteDomain<DDimY> const y_domain
             = y_domain_with_periodic_point.remove_last(ddc::DiscreteVector<DDimY>(1));
 
@@ -123,11 +123,8 @@ int main(int argc, char** argv)
     // Initialization of the global domain in time:
     // - the number of discrete time-points is equal to the number of
     //   steps + 1
-    ddc::DiscreteDomain<DDimT> const time_domain
-            = ddc::init_discrete_space<DDimT>(DDimT::init<DDimT>(
-                    ddc::Coordinate<T>(start_time),
-                    ddc::Coordinate<T>(end_time),
-                    nb_time_steps + 1));
+    ddc::DiscreteDomain<DDimT> const time_domain = ddc::create_uniform_point_sampling<
+            DDimT>(ddc::Coordinate<T>(start_time), ddc::Coordinate<T>(end_time), nb_time_steps + 1);
 
     ddc::DiscreteDomain<DDimX, DDimY> const xy_domain(x_domain, y_domain);
 
@@ -156,8 +153,8 @@ int main(int argc, char** argv)
     // time of the iteration where the last output happened
     ddc::DiscreteElement<DDimT> last_output = time_domain.front();
 
-    ddc::init_discrete_space<DDimFx>(ddc::init_fourier_space<DDimFx>(x_domain));
-    ddc::init_discrete_space<DDimFy>(ddc::init_fourier_space<DDimFy>(y_domain));
+    ddc::create_fourier_space<DDimFx>(x_domain);
+    ddc::create_fourier_space<DDimFy>(y_domain);
 
     ddc::DiscreteDomain<DDimFx, DDimFy> const k_mesh
             = ddc::fourier_mesh<DDimFx, DDimFy>(xy_domain, false);

diff --git a/examples/non_uniform_heat_equation.cpp b/examples/non_uniform_heat_equation.cpp
@@ -175,8 +175,8 @@ int main(int argc, char** argv)
 
     //! [build-domains]
     auto const [x_domain, ghosted_x_domain, x_pre_ghost, x_post_ghost]
-            = ddc::init_discrete_space<DDimX>(
-                    DDimX::init_ghosted<DDimX>(x_domain_vect, x_pre_ghost_vect, x_post_ghost_vect));
+            = ddc::create_non_uniform_point_sampling_ghosted<
+                    DDimX>(x_domain_vect, x_pre_ghost_vect, x_post_ghost_vect);
     //! [build-domains]
 
     ddc::DiscreteDomain<DDimX> const
@@ -195,8 +195,8 @@ int main(int argc, char** argv)
 
     //! [build-Y-domain]
     auto const [y_domain, ghosted_y_domain, y_pre_ghost, y_post_ghost]
-            = ddc::init_discrete_space<DDimY>(
-                    DDimY::init_ghosted<DDimY>(y_domain_vect, y_pre_ghost_vect, y_post_ghost_vect));
+            = ddc::create_non_uniform_point_sampling_ghosted<
+                    DDimY>(y_domain_vect, y_pre_ghost_vect, y_post_ghost_vect);
     //! [build-Y-domain]
 
     ddc::DiscreteDomain<DDimY> const
@@ -229,11 +229,8 @@ int main(int argc, char** argv)
     ddc::DiscreteVector<DDimT> const nb_time_steps(
             std::ceil((end_time - start_time) / max_dt) + .2);
 
-    ddc::DiscreteDomain<DDimT> const time_domain
-            = ddc::init_discrete_space<DDimT>(DDimT::init<DDimT>(
-                    ddc::Coordinate<T>(start_time),
-                    ddc::Coordinate<T>(end_time),
-                    nb_time_steps + 1));
+    ddc::DiscreteDomain<DDimT> const time_domain = ddc::create_uniform_point_sampling<
+            DDimT>(ddc::Coordinate<T>(start_time), ddc::Coordinate<T>(end_time), nb_time_steps + 1);
     //! [time-domain]
 
     //! [data allocation]

diff --git a/examples/uniform_heat_equation.cpp b/examples/uniform_heat_equation.cpp
@@ -96,11 +96,11 @@ int main(int argc, char** argv)
 
     //! [X-global-domain]
     auto const [x_domain, ghosted_x_domain, x_pre_ghost, x_post_ghost]
-            = ddc::init_discrete_space<DDimX>(DDimX::init_ghosted<DDimX>(
+            = ddc::create_uniform_point_sampling_ghosted<DDimX>(
                     ddc::Coordinate<X>(x_start),
                     ddc::Coordinate<X>(x_end),
                     ddc::DiscreteVector<DDimX>(nb_x_points),
-                    gwx));
+                    gwx);
     //! [X-global-domain]
 
     //! [X-domains]
@@ -114,11 +114,11 @@ int main(int argc, char** argv)
     ddc::DiscreteVector<DDimY> const gwy(1);
 
     auto const [y_domain, ghosted_y_domain, y_pre_ghost, y_post_ghost]
-            = ddc::init_discrete_space<DDimY>(DDimY::init_ghosted<DDimY>(
+            = ddc::create_uniform_point_sampling_ghosted<DDimY>(
                     ddc::Coordinate<Y>(y_start),
                     ddc::Coordinate<Y>(y_end),
                     ddc::DiscreteVector<DDimY>(nb_y_points),
-                    gwy));
+                    gwy);
 
     ddc::DiscreteDomain<DDimY> const
             y_post_mirror(y_post_ghost.front() - y_domain.extents(), y_post_ghost.extents());
@@ -139,11 +139,8 @@ int main(int argc, char** argv)
     //! [time-domain]
     ddc::DiscreteVector<DDimT> const nb_time_steps(std::ceil((end_time - start_time) / dt) + .2);
 
-    ddc::DiscreteDomain<DDimT> const time_domain
-            = ddc::init_discrete_space<DDimT>(DDimT::init<DDimT>(
-                    ddc::Coordinate<T>(start_time),
-                    ddc::Coordinate<T>(end_time),
-                    nb_time_steps + 1));
+    ddc::DiscreteDomain<DDimT> const time_domain = ddc::create_uniform_point_sampling<
+            DDimT>(ddc::Coordinate<T>(start_time), ddc::Coordinate<T>(end_time), nb_time_steps + 1);
     //! [time-domain]
 
     //! [data allocation]

diff --git a/include/ddc/discrete_space.hpp b/include/ddc/discrete_space.hpp
@@ -216,6 +216,16 @@ std::tuple<Arg0, Arg1, Args...> init_discrete_space(std::tuple<DDimImpl, Arg0, A
     return detail::extract_after(std::move(a), std::index_sequence_for<Arg0, Arg1, Args...>());
 }
 
+/** Initialize (emplace) a global singleton discrete space
+ *
+ * @param ddim_impl the impl discrete space object
+ */
+template <class DDim, class DDimImpl>
+void init_discrete_space_from_impl(DDimImpl&& ddim_impl)
+{
+    init_discrete_space<DDim>(std::forward<DDimImpl>(ddim_impl));
+}
+
 /**
  * @tparam DDim a discrete dimension
  * @return the discrete space instance associated with `DDim`.

diff --git a/include/ddc/experimental/single_discretization.hpp b/include/ddc/experimental/single_discretization.hpp
@@ -89,6 +89,14 @@ class SingleDiscretization : SingleDiscretizationBase
     };
 };
 
+template <class DDim>
+void create_single_discretization(Coordinate<typename DDim::continuous_dimension_type> origin)
+{
+    init_discrete_space_from_impl<DDim>(
+            typename SingleDiscretization<typename DDim::continuous_dimension_type>::
+                    template Impl<DDim, Kokkos::HostSpace>(origin));
+}
+
 template <class DDim>
 struct is_single_discretization : public std::is_base_of<SingleDiscretizationBase, DDim>::type
 {

diff --git a/include/ddc/kernels/fft.hpp b/include/ddc/kernels/fft.hpp
@@ -300,6 +300,47 @@ typename DDimFx::template Impl<DDimFx, Kokkos::HostSpace> init_fourier_space(
     return std::move(impl);
 }
 
+/**
+ * @brief Initialize a Fourier discrete dimension.
+ *
+ * Initialize the (1D) discrete space representing the Fourier discrete dimension associated
+ * to the (1D) mesh passed as argument. It is a N-periodic PeriodicSampling with a periodic window of width 2*pi/dx.
+ *
+ * This value comes from the Nyquist-Shannon theorem: the period of the spectral domain is N*dk = 2*pi/dx.
+ * Adding to this the relations dx = (xmax-xmin)/(N-1), and dk = (kmax-kmin)/(N-1), we get kmax-kmin = 2*pi*(N-1)^2/N/(xmax-xmin),
+ * which is used in the implementation (xmax, xmin, kmin and kmax are the centers of lower and upper cells inside a single period of the meshes).
+ *
+ * @tparam DDimFx A PeriodicSampling representing the Fourier discrete dimension.
+ * @tparam DDimX The type of the original discrete dimension.
+ *
+ * @param x_mesh The DiscreteDomain representing the (1D) original mesh.
+ *
+ * @see PeriodicSampling
+ */
+template <typename DDimFx, typename DDimX>
+void create_fourier_space(ddc::DiscreteDomain<DDimX> x_mesh)
+{
+    static_assert(
+            is_uniform_point_sampling_v<DDimX>,
+            "DDimX dimension must derive from UniformPointSampling");
+    static_assert(
+            is_periodic_sampling_v<DDimFx>,
+            "DDimFx dimension must derive from PeriodicSampling");
+    using CDimFx = typename DDimFx::continuous_dimension_type;
+    using CDimX = typename DDimX::continuous_dimension_type;
+    static_assert(
+            std::is_same_v<CDimFx, ddc::Fourier<CDimX>>,
+            "DDimX and DDimFx dimensions must be defined over the same continuous dimension");
+
+    DiscreteVectorElement const nx = get<DDimX>(x_mesh.extents());
+    double const lx = ddc::rlength(x_mesh);
+    create_periodic_sampling<DDimFx>(
+            ddc::Coordinate<CDimFx>(0),
+            ddc::Coordinate<CDimFx>(2 * (nx - 1) * (nx - 1) / (nx * lx) * Kokkos::numbers::pi),
+            ddc::DiscreteVector<DDimFx>(nx),
+            ddc::DiscreteVector<DDimFx>(nx));
+}
+
 /**
  * @brief Get the Fourier mesh.
  *

diff --git a/include/ddc/kernels/splines/bsplines_non_uniform.hpp b/include/ddc/kernels/splines/bsplines_non_uniform.hpp
@@ -385,6 +385,46 @@ class NonUniformBSplines : detail::NonUniformBSplinesBase
     };
 };
 
+/** @brief Constructs a non uniform spline basis by iterating over a range of break points from begin to end.
+ * The associated static attributes are available after this call.
+ *
+ * The provided break points describe the separation between the cells on which the polynomials
+ * comprising a spline are defined. They are used to build a set of knots. There are 2*degree more
+ * knots than break points. In the non-periodic case the knots are defined as follows:
+ * \f$ k_i = b_0 \forall 0 \leq i < d \f$
+ * \f$ k_{i+d} = b_i \forall 0 \leq i < n_b \f$
+ * \f$ k_{i+d+n_b} = b_{n_b-1} \forall 0 \leq i < d \f$
+ * where \f$d\f$ is the degree of the polynomials, and \f$n_b\f$ is the number of break points in the input pair of iterators. And in the periodic case:
+ * \f$ k_i = b_{n_b-1-d+i} \forall 0 \leq i < d \f$
+ * \f$ k_{i+d} = b_i \forall 0 \leq i \leq n_b \f$
+ * \f$ k_{i+d+n_b} = b_{i+1} \forall 0 \leq i < d \f$
+ *
+ * This constructor makes the knots accessible via a DiscreteSpace.
+ * @tparam DDim  The name of the discrete dimension
+ * @param breaks_begin The iterator which points at the beginning of the break points.
+ * @param breaks_end The iterator which points at the end of the break points.
+ */
+template <class DDim, class RandomIt>
+void create_non_uniform_bsplines(RandomIt const breaks_begin, RandomIt const breaks_end)
+{
+    ddc::init_discrete_space_from_impl<DDim>(
+            typename NonUniformBSplines<typename DDim::continuous_dimension_type, DDim::degree()>::
+                    template Impl<DDim, Kokkos::HostSpace>(breaks_begin, breaks_end));
+}
+
+/** @brief Constructs a non uniform spline basis using break points contained in a std::vector.
+ * The associated static attributes are available after this call.
+ *
+ * @tparam DDim  The name of the discrete dimension
+ * @param breaks The std::vector of the coordinates of break points.
+ */
+template <class DDim>
+void create_non_uniform_bsplines(
+        std::vector<ddc::Coordinate<typename DDim::continuous_dimension_type>> const& breaks)
+{
+    create_non_uniform_bsplines<DDim>(breaks.begin(), breaks.end());
+}
+
 template <class DDim>
 struct is_non_uniform_bsplines : public std::is_base_of<detail::NonUniformBSplinesBase, DDim>::type
 {
@@ -439,7 +479,7 @@ NonUniformBSplines<CDim, D>::Impl<DDim, MemorySpace>::Impl(
             knots[degree() + npoints() - 1 + i] = rmax;
         }
     }
-    ddc::init_discrete_space<knot_discrete_dimension_type>(knots);
+    ddc::create_non_uniform_point_sampling<knot_discrete_dimension_type>(knots);
 }
 
 template <class CDim, std::size_t D>