Skip to content

Commit

Permalink
Merge branch 'splines-linear-problem-3x3-blocks' into lapack-backend
Browse files Browse the repository at this point in the history
  • Loading branch information
blegouix committed Jun 26, 2024
2 parents 2b0ee29 + 8ab0cb2 commit 2597448
Show file tree
Hide file tree
Showing 7 changed files with 484 additions and 111 deletions.
4 changes: 2 additions & 2 deletions examples/non_uniform_heat_equation.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -17,8 +17,8 @@
//! [vector_generator]
std::vector<double> 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());
Expand Down
6 changes: 6 additions & 0 deletions include/ddc/kernels/splines/greville_interpolation_points.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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 <
Expand All @@ -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 <
Expand Down Expand Up @@ -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 <class Sampling>
Expand Down
4 changes: 4 additions & 0 deletions include/ddc/kernels/splines/knots_as_interpolation_points.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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 BSplines, ddc::BoundCond BcXmin, ddc::BoundCond BcXmax>
class KnotsAsInterpolationPoints
Expand Down
8 changes: 8 additions & 0 deletions include/ddc/kernels/splines/null_extrapolation_rule.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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 <class CoordType, class ChunkSpan>
KOKKOS_FUNCTION double operator()(CoordType, ChunkSpan) const
{
Expand Down
175 changes: 167 additions & 8 deletions include/ddc/kernels/splines/spline_evaluator.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -15,57 +15,93 @@

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,
class... IDimX>
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<evaluation_mesh_type>;

/// @brief The type of the whole domain representing evaluation points.
using batched_evaluation_domain_type = ddc::DiscreteDomain<IDimX...>;

/// @brief The type of the 1D spline domain corresponding to the dimension of interest.
using spline_domain_type = ddc::DiscreteDomain<bsplines_type>;

/**
* @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::type_seq_remove_t<
ddc::detail::TypeSeq<IDimX...>,
ddc::detail::TypeSeq<evaluation_mesh_type>>>;

/**
* @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::type_seq_replace_t<
ddc::detail::TypeSeq<IDimX...>,
ddc::detail::TypeSeq<evaluation_mesh_type>,
ddc::detail::TypeSeq<bsplines_type>>>;

/// @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;
Expand Down Expand Up @@ -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)
Expand All @@ -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 <class Layout, class... CoordsDims>
KOKKOS_FUNCTION double operator()(
ddc::Coordinate<CoordsDims...> const& coord_eval,
Expand All @@ -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 <class Layout1, class Layout2, class Layout3, class... CoordsDims>
void operator()(
ddc::ChunkSpan<double, batched_evaluation_domain_type, Layout1, memory_space> const
Expand Down Expand Up @@ -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 <class Layout, class... CoordsDims>
KOKKOS_FUNCTION double deriv(
ddc::Coordinate<CoordsDims...> const& coord_eval,
Expand All @@ -179,6 +307,24 @@ class SplineEvaluator
return eval_no_bc<eval_deriv_type>(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 <class Layout1, class Layout2, class Layout3, class... CoordsDims>
void deriv(
ddc::ChunkSpan<double, batched_evaluation_domain_type, Layout1, memory_space> const
Expand Down Expand Up @@ -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 <class Layout1, class Layout2>
void integrate(
ddc::ChunkSpan<double, batch_domain_type, Layout1, memory_space> const integrals,
Expand Down
Loading

0 comments on commit 2597448

Please sign in to comment.