Skip to content
GitLab
Commit c850b787 authored by COULAUD Olivier's avatar COULAUD Olivier
Browse files

Start the chebychev interpolator

parent c16a282e
Hide whitespace changes
Inline Side-by-side
experimental/examples/test_l2p.cpp
View file @ c850b787
......@@ -3,6 +3,7 @@
#include "scalfmm/container/particle_container.hpp"
#include "scalfmm/container/point.hpp"
#include "scalfmm/interpolation/uniform.hpp"
#include "scalfmm/interpolation/chebyshev.hpp"
#include "scalfmm/operators/fmmOperators.hpp"
#include "scalfmm/operators/generic/l2p.hpp"
#include "scalfmm/matrix_kernels/laplace.hpp"
......@@ -109,6 +110,7 @@ auto main([[maybe_unused]] int argc, [[maybe_unused]] char* argv[]) -> int
// Number of output p + force
static constexpr std::size_t nb_outputs_leaf{nb_outputs+dimension };
using interpolator_type = scalfmm::interpolation::uniform_interpolator<value_type, dimension, matrix_kernel_type>;
using interpolator_type2 = scalfmm::interpolation::chebyshev_interpolator<value_type, dimension, matrix_kernel_type>;
using far_field_type = scalfmm::operators::far_field_operator<interpolator_type, true>;
// pos, charge, pot + forces no variables
......@@ -189,6 +191,13 @@ auto main([[maybe_unused]] int argc, [[maybe_unused]] char* argv[]) -> int
static_cast<std::size_t>(tree_height),
box.width(0));
far_field_type far_uniform(uniform);
std::cout << " Start Cheb " << std::endl;
interpolator_type2 cheb(matrix_kernel_type{}, order,
static_cast<std::size_t>(tree_height),
box.width(0));
auto roots2_1d = cheb.roots(static_cast<std::size_t>(order));
std::cout << "roots: " << roots2_1d << std::endl;
std::exit(EXIT_FAILURE);
///
/////////////////////////////////////////////////////////////////////////////
......@@ -202,6 +211,9 @@ auto main([[maybe_unused]] int argc, [[maybe_unused]] char* argv[]) -> int
//
auto roots_1d = uniform.roots(static_cast<std::size_t>(order));
std::cout << "roots: " << roots_1d <<std::endl;
auto roots = scalfmm::tensor::generate_meshgrid<dimension>(roots_1d);
const auto roots_x = xt::flatten(std::get<0>(roots));
const auto roots_y = xt::flatten(std::get<1>(roots));
......@@ -209,7 +221,6 @@ auto main([[maybe_unused]] int argc, [[maybe_unused]] char* argv[]) -> int
std::cout << "Roots_x: " << roots_x << std::endl;
std::cout << "Roots_y: " << roots_y << std::endl;
std::cout << "Roots_z: " << roots_z << std::endl;
// std::cout << "roots: " << roots_1d <<std::endl;
int pos = 0, global_idx = 0;
scalfmm::component::for_each_cell_leaf(
......
experimental/include/scalfmm/interpolation/chebyshev.hpp
View file @ c850b787
// --------------------------------
// --------------------------------
// See LICENCE file at project root
// File : interpolation/chebyshev_interpolator.hpp
// File : interpolation/chebyshev.hpp
// --------------------------------
#ifndef SCALFMM_INTERPOLATION_CHEBYSHEV_INTERPOLATOR_HPP
#define SCALFMM_INTERPOLATION_CHEBYSHEV_INTERPOLATOR_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/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 <cmath>
#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 <scalfmm/interpolation/interpolator.hpp>
#include <scalfmm/simd/utils.hpp>
#include <xsimd/math/xsimd_math.hpp>
#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
{
template<typename T>
struct chebyshev_roots : roots_callable<T>
template<typename ValueType, std::size_t Dimension, typename FarFieldMatrixKernel>
struct chebyshev_interpolator
: public interpolator<chebyshev_interpolator<ValueType, Dimension, FarFieldMatrixKernel>>
{
using array_type = typename roots_callable<T>::array_type;
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 complex_type = std::complex<value_type>;
using k_tensor_type = xt::xarray<complex_type>;
using interaction_matrix_type = xt::xtensor_fixed<k_tensor_type, xt::xshape<kn, km>>;
using base_type::base_type;
using buffer_inner_type = xt::xarray<std::complex<value_type>>;
using buffer_shape_type = xt::xshape<kn>;
using buffer_type = xt::xtensor_fixed<buffer_inner_type, buffer_shape_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;
inline array_type operator()(std::size_t order) override
[[nodiscard]] inline auto roots_impl(std::size_t order) const
{
array_type roots = xt::arange(T(1.0), T(order + 1));
const value_type coeff = value_type(3.14159265358979323846264338327950) / (value_type(order));
// xfunction xt::cos
std::cout << " cos(n) " << (xt::cos(coeff * (xt::arange(int(order - 1), int(-1), -1) + 0.5))) << std::endl;
// return xt::cos(coeff * (xt::arange(int(order - 1), int(-1), -1) + 0.5));
return xt::linspace(value_type(-1.), value_type(1), order);
}
template<typename ComputationType>
[[nodiscard]] inline auto polynomials_impl(ComputationType x, std::size_t order, std::size_t n) const
{
using value_type = ComputationType;
// assert(xsimd::any(xsimd::abs(x) - 1. < 10. * std::numeric_limits<value_type>::epsilon()));
const value_type two_const{2.0};
simd::compare(x);
auto cheb_formula = [order](auto& e) {
e = xsimd::cos(
((T(2.) * e - T(1.)) / (T(2.) * T(order))) *
T(3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651e+00));
};
// Specific precomputation of scale factor
// in order to minimize round-off errors
// NB: scale factor could be hardcoded (just as the roots)
value_type scale{};
const int omn{int(order) - int(n) - 1};
std::for_each(roots.begin(), roots.end(), cheb_formula);
if(omn % 2 != 0)
{
scale = value_type(-1.); // (-1)^(n-1-(k+1)+1)=(-1)^(omn-1)
}
else
{
scale = value_type(1.);
}
scale /= xsimd::pow(two_const, int(order) - 1) * math::factorial<value_type>(int(n)) *
math::factorial<value_type>(omn);
// compute L
value_type L{1.};
const value_type coeff = value_type(int(order) - 1);
for(std::size_t m = 0; m < order; ++m)
{
if(m != n)
{
// previous version with risks of round-off error
// L *= (x-FUnifRoots<order>::roots[m])/(FUnifRoots<order>::roots[n]-FUnifRoots<order>::roots[m]);
return roots;
// new version (reducing round-off)
// regular grid on [-1,1] (h simplifies, only the size of the domain and a remains i.e. 2. and -1.)
L *= ((value_type(int(order) - 1) * (x + value_type(1.))) - (two_const * value_type(int(m))));
// L *= ((coeff * (x + value_type(1.))) - (two_const * value_type(int(m))));
}
}
L *= scale;
return L;
}
inline T operator()(std::size_t order, std::size_t n)
template<typename ComputationType>
[[nodiscard]] inline auto derivative_impl(ComputationType x, std::size_t order, std::size_t n) const
{
return xsimd::cos(
((T(2.) * T(n) - T(1.)) / (T(2.) * T(order))) *
T(3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651e+00));
using value_type = ComputationType;
// assert(xsimd::any(xsimd::abs(x) - 1. < 10. * std::numeric_limits<value_type>::epsilon()));
const value_type two_const{2.0};
simd::compare(x);
// optimized variant
value_type NdL(0.); // init numerator
value_type DdL(1.); // init denominator
value_type tmpNdL{};
auto roots(roots_impl(order));
for(unsigned int p = 0; p < order; ++p)
{
if(p != n)
{
tmpNdL = value_type(1.);
for(unsigned int m = 0; m < order; ++m)
{
if(m != n && m != p)
{
tmpNdL *= x - roots.at(m);
}
}
NdL += tmpNdL;
DdL *= roots.at(n) - roots.at(p);
} // endif
} // p
return NdL / DdL;
}
};
template<typename T>
struct chebyshev_polynomials : polynomials_callable<T>
{
using value_type = T;
// explicit chebyshev_polynomials(std::size_t order)
//{
// m_roots(m_roots_function(order));
// m_t_of_roots.resize({order, order});
// for(std::size_t o = 1; o < order; ++o)
// {
// for(std::size_t oo = 1; oo < order; ++oo)
// {
// m_t_of_roots(o, oo) = p_of_first_kind(m_roots[oo], o);
// }
// }
//}
inline value_type p_of_first_kind(value_type x, std::size_t n)
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
{
simd::compare(x);
return value_type(xsimd::cos(n * xsimd::acos(x)));
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);
}
inline value_type operator()(value_type x, const std::size_t order, const std::size_t n) override
[[nodiscard]] inline auto buffer_initialization_impl() const -> buffer_type
{
value_type s{value_type(1.) / order};
std::vector<std::size_t> shape(dimension, 2 * this->order() - 1);
shape.at(dimension - 1) = this->order();
return buffer_type(buffer_shape_type{}, buffer_inner_type(shape, 0.));
}
for(std::size_t o = 1; o < order; ++o)
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)
{
auto t_of_x = p_of_first_kind(x, o);
auto t_of_roots = p_of_first_kind(m_roots_function(order, n), o);
s += value_type(2.) / order * t_of_x * t_of_roots;
mtilde.at(m).resize(shape_nnodes);
}
return s;
}
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.));
});
private:
xt::xarray<value_type> m_roots{};
xt::xarray<value_type> m_t_of_roots{};
chebyshev_roots<value_type> m_roots_function{};
};
auto views_in_padded = meta::apply(
padded_tensors, [&order](auto& t) { return tensor::get_view<dimension>(t, xt::range(0, order)); });
template<typename ValueType, std::size_t Dimension>
struct chebyshev_interpolator : public interpolator<ValueType, Dimension>
views_in_padded = meta::apply(multipoles, [](auto const& r) { return r; });
{
public:
using base_type = interpolator<ValueType, Dimension>;
using size_type = typename base_type::size_type;
using value_type = typename base_type::value_type;
using base_type::base_type;
meta::repeat_indexed(
[&mtilde](std::size_t m, auto const& r) { mtilde.at(m) = xt::fftw::rfftn<dimension>(r); },
padded_tensors);
}
template<typename T>
using roots_f_type = chebyshev_roots<T>;
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();
template<typename T>
using polynomials_f_type = chebyshev_polynomials<T>;
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);
}
}
explicit chebyshev_interpolator(size_type order, size_type tree_height = 3,
value_type root_cell_width = value_type(1.),
value_type cell_width_extension = value_type(0.))
: base_type(order, chebyshev_roots<value_type>{}, chebyshev_polynomials<value_type>{}, tree_height,
root_cell_width, cell_width_extension)
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);
}
};
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<std::complex<value_type>>;
using buffer_inner_type = xt::xarray<std::complex<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_type = xt::xtensor_fixed<buffer_inner_type, xt::xshape<far_field_matrix_kernel_type::kn>>;
};
} // namespace scalfmm::interpolation
#endif // SCALFMM_INTERPOLATION_CHEBYSHEV_INTERPOLATOR_HPP
#endif // SCALFMM_INTERPOLATION_CHEBYSHEV_HPP
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment