Skip to content

GitLab

Commit bed26e0a authored by ESTERIE Pierre's avatar ESTERIE Pierre
Browse files

Fix Chebychev initialisation

parent a17c2d6f
Hide whitespace changes
Inline Side-by-side
experimental/include/scalfmm/interpolation/chebyshev.hpp
View file @ bed26e0a
// --------------------------------
// --------------------------------
// See LICENCE file at project root
// File : interpolation/chebyshev.hpp
// --------------------------------
#ifndef SCALFMM_INTERPOLATION_CHEBYSHEV_HPP
#define SCALFMM_INTERPOLATION_CHEBYSHEV_HPP
#include "scalfmm/container/point.hpp"
#include "scalfmm/container/variadic_adaptor.hpp"
#include "scalfmm/interpolation/generate_circulent.hpp"
#include "scalfmm/interpolation/grid_storage.hpp"
#include "scalfmm/interpolation/interpolator.hpp"
#include "scalfmm/interpolation/mapping.hpp"
#include "scalfmm/interpolation/traits.hpp"
#include "scalfmm/meta/traits.hpp"
#include "scalfmm/meta/utils.hpp"
#include "scalfmm/simd/memory.hpp"
#include "scalfmm/simd/utils.hpp"
#include "scalfmm/tools/colorized.hpp"
#include "scalfmm/utils/fftw.hpp"
#include "scalfmm/utils/math.hpp"
#include "scalfmm/utils/tensor.hpp"
#include <algorithm>
#include <array>
#include <bits/c++config.h>
#include <cassert>
#include <cmath>
#include <complex>
#include <iterator>
#include <limits>
#include <memory>
#include <tuple>
#include <type_traits>
#include <utility>
#include <xsimd/memory/xsimd_load_store.hpp>
#include <xsimd/xsimd.hpp>
#include <xtensor-fftw/basic.hpp>
#include <xtensor/xarray.hpp>
#include <xtensor/xbuilder.hpp>
#include <xtensor/xcomplex.hpp>
#include <xtensor/xio.hpp>
#include <xtensor/xmanipulation.hpp>
#include <xtensor/xmath.hpp>
#include <xtensor/xnoalias.hpp>
#include <xtensor/xslice.hpp>
#include <xtensor/xtensor.hpp>
#include <xtensor/xtensor_forward.hpp>
#include <xtensor/xtensor_simd.hpp>
#include <xtensor/xvectorize.hpp>
#include <xtl/xclosure.hpp>
#include <xtl/xcomplex.hpp>
namespace scalfmm::interpolation
{
struct empty
{
};
template<typename ValueType, std::size_t Dimension, typename FarFieldMatrixKernel>
struct chebyshev_interpolator
: public interpolator<chebyshev_interpolator<ValueType, Dimension, FarFieldMatrixKernel>>
{
public:
using value_type = ValueType;
static constexpr std::size_t dimension = Dimension;
using matrix_kernel_type = FarFieldMatrixKernel;
static constexpr std::size_t kn = matrix_kernel_type::kn;
static constexpr std::size_t km = matrix_kernel_type::km;
using base_type = interpolator<chebyshev_interpolator<value_type, dimension, matrix_kernel_type>>;
using k_tensor_type = xt::xarray<value_type>;
using storage_type = component::grid_storage<value_type, dimension, km, kn>;
using interaction_matrix_type = xt::xtensor_fixed<k_tensor_type, xt::xshape<kn, km>>;
using base_type::base_type;
using buffer_inner_type = empty;
// using buffer_shape_type = xt::xshape<kn>;
using buffer_type = empty;
chebyshev_interpolator() = default;
chebyshev_interpolator(chebyshev_interpolator const&) = default;
chebyshev_interpolator(chebyshev_interpolator&&) noexcept = default;
auto operator=(chebyshev_interpolator const&) -> chebyshev_interpolator& = default;
auto operator=(chebyshev_interpolator&&) noexcept -> chebyshev_interpolator& = default;
~chebyshev_interpolator() = default;
///
/// \brief roots_impl Chebyshev roots in [-1,1]
///
/// computed as \f$\bar x_n = \cos\left(\frac{\pi}{2}\frac{2n-1}{\ell}\right)\f$ for \f$n=1,\dots,\ell\f$
/// \param order
/// \return the Chebyshev roots
///
[[nodiscard]] inline auto roots_impl(std::size_t order) const -> xt::xarray<value_type>
{
const value_type coeff = value_type(3.14159265358979323846264338327950) / (value_type(order));
return xt::cos(coeff * (xt::arange(int(order - 1), int(-1), -1) + 0.5));
}
/**
* S_k(x) Lagrange fonction based on Chebyshev polynomials of first kind \f$
*
*for \f$p = order+1 \f$ the interpolator S assocoated to point n
* \f$ S_n(x) = \frac{1}{p} + \frac{2}{p} \sum_{k=1}{p-1}{ T_k(x)T_k(x_n)}\f$
* then we have \f$ S_n(x_m) =\delta_{n,m} \f$
*\todo expensive method !!
*
* @param[in] n index
* @param[in] x coordinate in [-1,1]
* @return function value
*/
template<typename ComputationType>
[[nodiscard]] inline auto polynomials_impl(ComputationType x, std::size_t n) const
{
using value_type = ComputationType;
simd::compare(x);
auto order{this->order()};
//
// value_type* T_of_x = new value_type[order];
auto coeff{value_type(2.0) / value_type(order)};
value_type L{0.5};
for(unsigned int o = 1; o < order; ++o)
{
L += T(o, x) * m_T_of_roots[o * order + n];
}
return coeff * L;
}
/**
* d/dx S_n(x) gradient of the Lagrange fonction based on Chebyshev
*polynomials of first kind \f$
*
*for \f$p = order+1 \f$ the interpolator S assocoated to point n
* \f$ d/dx S_n(x) = \frac{2}{p} \sum_{k=1}{p-1}{ T_k(x_n) d/dx T_k(x)}\f$
* then
* \f$ d/dx S_n(x) = \frac{2}{p} \sum_{k=1}{p}{ T_k(x_n) k U_{k-1}(x)}\f$
*
* \todo expensive method !!
*
* @param[in] n index
* @param[in] x coordinate in [-1,1]
* @return function value
*/
template<typename ComputationType>
[[nodiscard]] inline auto derivative_impl(ComputationType x, std::size_t n) const
{
using value_type = ComputationType;
simd::compare(x);
auto order{this->order()};
//
auto coeff{value_type(2.0) / value_type(order)};
// U(0, x) = 1
value_type L{m_T_of_roots[order + n]};
//{value_type(0)};
for(unsigned int k = 2; k < order; ++k)
{
L += k * U(k - 1, x) * m_T_of_roots[k * order + n];
}
return coeff * L;
}
inline auto generate_interactions_matrixes_impl(value_type width, std::size_t tree_height) -> void {}
template<typename TensorViewX, typename TensorViewY>
[[nodiscard]] inline auto generate_matrix_k_impl(TensorViewX&& X, TensorViewY&& Y, std::size_t order,
matrix_kernel_type const& far_field, std::size_t n,
std::size_t m) const -> k_tensor_type
{
build<value_type, dimension, dimension> build_c{};
auto C = build_c(std::forward<TensorViewX>(X), std::forward<TensorViewY>(Y), order, far_field, n, m);
return xt::fftw::rfftn<dimension>(C);
}
[[nodiscard]] inline auto buffer_initialization_impl() const -> buffer_type { return empty{}; }
template<typename Cell>
void apply_multipoles_preprocessing_impl(Cell& current_cell, std::size_t order) const
{
// using multipoles_container = typename Cell::multipoles_container;
// using fft_type = container::get_variadic_adaptor_t<xt::xarray<typename
// multipoles_container::value_type>,
// Cell::inputs_size>;
// auto const& multipoles = current_cell.cmultipoles();
// auto& mtilde = current_cell.transformed_multipoles();
// std::vector shape_nnodes(Cell::dimension, order);
// std::vector shape_transformed(Cell::dimension, std::size_t(2 * order - 1));
// // Resize buffers to store fft results
// for(std::size_t m = 0; m < km; ++m)
// {
// mtilde.at(m).resize(shape_nnodes);
// }
// fft_type padded_tensors;
// meta::for_each(padded_tensors, [&shape_transformed](auto& t) { t.resize(shape_transformed);
// }); meta::for_each(padded_tensors, [](auto& t) {
// std::fill(t.begin(), t.end(), typename multipoles_container::value_type(0.));
// });
// auto views_in_padded = meta::apply(
// padded_tensors, [&order](auto& t) { return tensor::get_view<dimension>(t, xt::range(0,
// order)); });
// views_in_padded = meta::apply(multipoles, [](auto const& r) { return r; });
// meta::repeat_indexed(
// [&mtilde](std::size_t m, auto const& r) { mtilde.at(m) = xt::fftw::rfftn<dimension>(r); },
// padded_tensors);
}
template<typename Cell, typename TupleScaleFactor>
void apply_m2l_impl(Cell const& source_cell, Cell& target_cell, interaction_matrix_type const& k,
TupleScaleFactor scale_factor, std::size_t order, std::size_t tree_level,
[[maybe_unused]] buffer_type& products) const
{
// auto const& transformed_multipoles = source_cell.ctransformed_multipoles();
// for(std::size_t m = 0; m < km; ++m)
// {
// meta::repeat_indexed(
// [&k, &products, &transformed_multipoles, m](std::size_t n, auto const& s_f) {
// products.at(n) += transformed_multipoles.at(m) * k.at(n, m) * s_f;
// },
// scale_factor);
// }
}
template<typename Cell>
void apply_multipoles_postprocessing_impl(Cell& current_cell,
[[maybe_unused]] buffer_type const& products) const
{
// auto order{current_cell.order()};
// auto& target_expansion = current_cell.locals();
// meta::repeat_indexed(
// [&products, order](std::size_t n, auto& expansion) {
// auto view_for_update =
// xt::eval(tensor::get_view<dimension>(
// xt::fftw::irfftn<dimension>(products.at(n),
// true), xt::range(0, order)));
// expansion += view_for_update;
// },
// target_expansion);
}
inline void set_private_members_impl()
{
auto order = this->order();
auto roots = this->roots(this->order());
// initialize chebyshev polynomials of root nodes: T_o(x_j)
/// m_T_of_roots[o,j] = T_o(x_j)
m_T_of_roots = new value_type[order * order];
std::cout << "size: " << order * order << std::endl;
for(unsigned int o = 0; o < order; ++o)
{
for(unsigned int j = 0; j < order; ++j)
{
m_T_of_roots[o * order + j] = T(o, roots[j]);
}
}
std::cout << " End of set_private_members_impl " << std::endl;
}
private:
///
/// \brief Chebyshev polynomials of first kind \f$ T_n(x) = \cos(n \arccos(x)) \f$
/// * \arccos(x))}{\sqrt{1-x^2}} \f$ are needed.
/// \param[in] n index
/// \param[in] x coordinate in [-1,1]
/// \return function value U_n(x)
template<typename ComputationType>
static ComputationType T(const unsigned int n, ComputationType x)
{
simd::compare(x);
return ComputationType(cos(n * acos(x)));
}
///
/// \brief U the derivation of the Chebyshev polynomials of first kind
///
/// For the derivation of the Chebyshev polynomials of first kind
/// \f$ \frac{\mathrm{d} T_n(x)}{\mathrm{d}x} = n U_{n-1}(x) \f$ the Chebyshev
/// polynomials of second kind \f$ U_{n-1}(x) = \frac{\sin(n \arccos(x))}{\sqrt{1-x^2}} \f$
/// are needed.
/// \param[in] n index
/// \param[in] x coordinate in [-1,1]
/// \return function value U_n(x)
///
template<typename ComputationType>
static ComputationType U(const unsigned int n, ComputationType x)
{
simd::compare(x);
return ComputationType((sin((n + 1) * acos(x))) / sqrt(1 - x * x));
}
/// Store Tn(x_i)
ValueType* m_T_of_roots; ///< The chebyshev polynomials of root nodes: m_T_of_roots[o,j] = T_o(x_j)
};
template<typename ValueType, std::size_t Dimension, typename FarFieldMatrixKernel>
struct interpolator_traits<chebyshev_interpolator<ValueType, Dimension, FarFieldMatrixKernel>>
{
using value_type = ValueType;
static constexpr std::size_t dimension = Dimension;
using far_field_matrix_kernel_type = FarFieldMatrixKernel;
using base_type = interpolator<chebyshev_interpolator<value_type, dimension, far_field_matrix_kernel_type>>;
using k_tensor_type = xt::xarray<value_type>;
using interaction_matrix_type =
xt::xtensor_fixed<k_tensor_type,
xt::xshape<far_field_matrix_kernel_type::kn, far_field_matrix_kernel_type::km>>;
using buffer_inner_type = empty;
using buffer_type = empty;
};
} // namespace scalfmm::interpolation
#include "scalfmm/interpolation/chebyshev/chebyshev_interpolator.hpp"
#endif // SCALFMM_INTERPOLATION_CHEBYSHEV_HPP
#endif // SCALFMM_INTERPOLATION_UNIFORM_HPP
experimental/include/scalfmm/interpolation/chebyshev/chebyshev_interpolator.hpp 0 → 100644
View file @ bed26e0a
// --------------------------------
// See LICENCE file at project root
// File : interpolation/chebyshev.hpp
// --------------------------------
#ifndef SCALFMM_INTERPOLATION_CHEBYSHEV_CHEBYSHEV_INTERPOLATOR_HPP
#define SCALFMM_INTERPOLATION_CHEBYSHEV_CHEBYSHEV_INTERPOLATOR_HPP
#include "scalfmm/container/point.hpp"
#include "scalfmm/container/variadic_adaptor.hpp"
#include "scalfmm/interpolation/generate_circulent.hpp"
#include "scalfmm/interpolation/grid_storage.hpp"
#include "scalfmm/interpolation/interpolator.hpp"
#include "scalfmm/interpolation/mapping.hpp"
#include "scalfmm/interpolation/traits.hpp"
#include "scalfmm/meta/traits.hpp"
#include "scalfmm/meta/utils.hpp"
#include "scalfmm/simd/memory.hpp"
#include "scalfmm/simd/utils.hpp"
#include "scalfmm/tools/colorized.hpp"
#include "scalfmm/utils/fftw.hpp"
#include "scalfmm/utils/math.hpp"
#include "scalfmm/utils/tensor.hpp"
#include <algorithm>
#include <array>
#include <bits/c++config.h>
#include <cassert>
#include <cmath>
#include <complex>
#include <iterator>
#include <limits>
#include <memory>
#include <tuple>
#include <type_traits>
#include <utility>
#include <xsimd/memory/xsimd_load_store.hpp>
#include <xsimd/xsimd.hpp>
#include <xtensor-fftw/basic.hpp>
#include <xtensor/xarray.hpp>
#include <xtensor/xbuilder.hpp>
#include <xtensor/xcomplex.hpp>
#include <xtensor/xio.hpp>
#include <xtensor/xmanipulation.hpp>
#include <xtensor/xmath.hpp>
#include <xtensor/xnoalias.hpp>
#include <xtensor/xslice.hpp>
#include <xtensor/xtensor.hpp>
#include <xtensor/xtensor_forward.hpp>
#include <xtensor/xtensor_simd.hpp>
#include <xtensor/xvectorize.hpp>
#include <xtl/xclosure.hpp>
#include <xtl/xcomplex.hpp>
namespace scalfmm::interpolation
{
struct empty
{
};
template<typename ValueType, std::size_t Dimension, typename FarFieldMatrixKernel>
struct chebyshev_interpolator
: public interpolator<chebyshev_interpolator<ValueType, Dimension, FarFieldMatrixKernel>>
{
public:
using value_type = ValueType;
static constexpr std::size_t dimension = Dimension;
using matrix_kernel_type = FarFieldMatrixKernel;
using size_type = std::size_t;
static constexpr std::size_t kn = matrix_kernel_type::kn;
static constexpr std::size_t km = matrix_kernel_type::km;
using base_type = interpolator<chebyshev_interpolator<value_type, dimension, matrix_kernel_type>>;
using k_tensor_type = xt::xarray<value_type>;
using storage_type = component::grid_storage<value_type, dimension, km, kn>;
using interaction_matrix_type = xt::xtensor_fixed<k_tensor_type, xt::xshape<kn, km>>;
using buffer_inner_type = empty;
using buffer_type = empty;
using base_type::base_type;
chebyshev_interpolator() = default;
chebyshev_interpolator(chebyshev_interpolator const&) = default;
chebyshev_interpolator(chebyshev_interpolator&&) noexcept = default;
auto operator=(chebyshev_interpolator const&) -> chebyshev_interpolator& = default;
auto operator=(chebyshev_interpolator&&) noexcept -> chebyshev_interpolator& = default;
~chebyshev_interpolator() = default;
chebyshev_interpolator(matrix_kernel_type&& far_field, size_type order, size_type tree_height = 3,
value_type root_cell_width = value_type(1.),
value_type cell_width_extension = value_type(0.), std::size_t separation_criterion = 1)
: base_type(std::forward<matrix_kernel_type>(far_field), order, tree_height, root_cell_width,
cell_width_extension, separation_criterion, true)
{
set_m_T_of_roots(order);
base_type::initialize(root_cell_width, tree_height);
}
chebyshev_interpolator(matrix_kernel_type const& far_field, size_type order, size_type tree_height = 3,
value_type root_cell_width = value_type(1.),
value_type cell_width_extension = value_type(0.), std::size_t separation_criterion = 1)
: base_type(std::forward<matrix_kernel_type>(far_field), order, tree_height, root_cell_width,
cell_width_extension, separation_criterion, true)
{
set_m_T_of_roots(order);
base_type::initialize(root_cell_width, tree_height);
}
///
/// \brief roots_impl Chebyshev roots in [-1,1]
///
/// computed as \f$\bar x_n = \cos\left(\frac{\pi}{2}\frac{2n-1}{\ell}\right)\f$ for \f$n=1,\dots,\ell\f$
/// \param order
/// \return the Chebyshev roots
///
[[nodiscard]] inline auto roots_impl(std::size_t order) const -> xt::xarray<value_type>
{
const value_type coeff = value_type(3.14159265358979323846264338327950) / (value_type(order));
return xt::cos(coeff * (xt::arange(int(order - 1), int(-1), -1) + 0.5));
}
/**
* S_k(x) Lagrange fonction based on Chebyshev polynomials of first kind \f$
*
*for \f$p = order+1 \f$ the interpolator S assocoated to point n
* \f$ S_n(x) = \frac{1}{p} + \frac{2}{p} \sum_{k=1}{p-1}{ T_k(x)T_k(x_n)}\f$
* then we have \f$ S_n(x_m) =\delta_{n,m} \f$
*\todo expensive method !!
*
* @param[in] n index
* @param[in] x coordinate in [-1,1]
* @return function value
*/
template<typename ComputationType>
[[nodiscard]] inline auto polynomials_impl(ComputationType x, std::size_t n) const
{
using value_type = ComputationType;
simd::compare(x);
auto order{this->order()};
auto coeff{value_type(2.0) / value_type(order)};
value_type L{0.5};
for(unsigned int o = 1; o < order; ++o)
{
L += T(o, x) * m_T_of_roots.at(o, n);
}
return coeff * L;
}
/**
* d/dx S_n(x) gradient of the Lagrange fonction based on Chebyshev
*polynomials of first kind \f$
*
*for \f$p = order+1 \f$ the interpolator S assocoated to point n
* \f$ d/dx S_n(x) = \frac{2}{p} \sum_{k=1}{p-1}{ T_k(x_n) d/dx T_k(x)}\f$
* then
* \f$ d/dx S_n(x) = \frac{2}{p} \sum_{k=1}{p}{ T_k(x_n) k U_{k-1}(x)}\f$
*
* \todo expensive method !!
*
* @param[in] n index
* @param[in] x coordinate in [-1,1]
* @return function value
*/
template<typename ComputationType>
[[nodiscard]] inline auto derivative_impl(ComputationType x, std::size_t n) const
{
using value_type = ComputationType;
simd::compare(x);
auto order{this->order()};
//
auto coeff{value_type(2.0) / value_type(order)};
// U(0, x) = 1
value_type L{m_T_of_roots.at(1, n)};
//{value_type(0)};
for(unsigned int k = 2; k < order; ++k)
{
L += k * U(k - 1, x) * m_T_of_roots.at(k, n);
}
return coeff * L;
}
inline auto generate_interactions_matrixes_impl(value_type width, std::size_t tree_height) -> void {}
template<typename TensorViewX, typename TensorViewY>
[[nodiscard]] inline auto generate_matrix_k_impl(TensorViewX&& X, TensorViewY&& Y, std::size_t order,
matrix_kernel_type const& far_field, std::size_t n,
std::size_t m) const -> k_tensor_type
{
build<value_type, dimension, dimension> build_c{};
auto C = build_c(std::forward<TensorViewX>(X), std::forward<TensorViewY>(Y), order, far_field, n, m);
return xt::fftw::rfftn<dimension>(C);
}
[[nodiscard]] inline auto buffer_initialization_impl() const -> buffer_type { return empty{}; }
template<typename Cell>
inline auto apply_multipoles_preprocessing_impl(Cell& current_cell, std::size_t order) const -> void
{
// using multipoles_container = typename Cell::multipoles_container;
// using fft_type = container::get_variadic_adaptor_t<xt::xarray<typename
// multipoles_container::value_type>,
// Cell::inputs_size>;
// auto const& multipoles = current_cell.cmultipoles();
// auto& mtilde = current_cell.transformed_multipoles();
// std::vector shape_nnodes(Cell::dimension, order);
// std::vector shape_transformed(Cell::dimension, std::size_t(2 * order - 1));
// // Resize buffers to store fft results
// for(std::size_t m = 0; m < km; ++m)
// {
// mtilde.at(m).resize(shape_nnodes);
// }
// fft_type padded_tensors;
// meta::for_each(padded_tensors, [&shape_transformed](auto& t) { t.resize(shape_transformed);
// }); meta::for_each(padded_tensors, [](auto& t) {
// std::fill(t.begin(), t.end(), typename multipoles_container::value_type(0.));
// });
// auto views_in_padded = meta::apply(
// padded_tensors, [&order](auto& t) { return tensor::get_view<dimension>(t, xt::range(0,
// order)); });
// views_in_padded = meta::apply(multipoles, [](auto const& r) { return r; });
// meta::repeat_indexed(
// [&mtilde](std::size_t m, auto const& r) { mtilde.at(m) = xt::fftw::rfftn<dimension>(r); },
// padded_tensors);
}
template<typename Cell, typename TupleScaleFactor>
inline auto apply_m2l_impl(Cell const& source_cell, Cell& target_cell, interaction_matrix_type const& k,
TupleScaleFactor scale_factor, std::size_t order, std::size_t tree_level,
[[maybe_unused]] buffer_type& products) const -> void
{
// auto const& transformed_multipoles = source_cell.ctransformed_multipoles();
// for(std::size_t m = 0; m < km; ++m)
// {
// meta::repeat_indexed(
// [&k, &products, &transformed_multipoles, m](std::size_t n, auto const& s_f) {
// products.at(n) += transformed_multipoles.at(m) * k.at(n, m) * s_f;
// },
// scale_factor);
// }
}
template<typename Cell>
inline auto apply_multipoles_postprocessing_impl(Cell& current_cell,
[[maybe_unused]] buffer_type const& products) const -> void
{
// auto order{current_cell.order()};
// auto& target_expansion = current_cell.locals();
// meta::repeat_indexed(
// [&products, order](std::size_t n, auto& expansion) {
// auto view_for_update =
// xt::eval(tensor::get_view<dimension>(