Commit 3f80e77c authored by COULAUD Olivier's avatar COULAUD Olivier
Browse files

Chebychev L2P works

parent aacab1f2
......@@ -163,14 +163,13 @@ auto run(const std::string& title, const std::string& input_file, const int &tre
auto roots_1d = interpolator.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));
const auto roots_z = xt::flatten(std::get<2>(roots));
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_x: " << roots_x << std::endl;
// std::cout << "Roots_y: " << roots_y << std::endl;
// std::cout << "Roots_z: " << roots_z << std::endl;
int pos = 0, global_idx = 0;
scalfmm::component::for_each_cell_leaf(
......@@ -201,19 +200,22 @@ auto run(const std::string& title, const std::string& input_file, const int &tre
// std::cout << " pos_z: " << pos_z << std::endl;
//
// Set the local expansion by f(pos_x, pos_y, pos_z)
auto& local_expansion = std::get<0>(cell.locals());
auto flatten_view = xt::flatten(local_expansion);
for(int i = 0; i < local_expansion.size(); ++i)
{
flatten_view.at(i) = func(pos_x[i], pos_y[i], pos_z[i]);
}
//
// Apply the l2p operator
/// farfield operator
///
scalfmm::operators::apply_l2p_der<far_field_type::compute_gradient>(far_field, cell.clocals(), leaf, order);
#ifndef toto
// for the next call nb_outputs_leaf = nb_outputs + dimension;
//scalfmm::operators::apply_l2p(far_uniform, cell.clocals(), leaf, order);
//scalfmm::operators::apply_l2p(far_field, cell.clocals(), leaf, order);
//
// Check the value on the particles
......@@ -253,9 +255,10 @@ auto run(const std::string& title, const std::string& input_file, const int &tre
}
//
#endif
++pos;
});
}
return 0;
}
......
......@@ -53,6 +53,9 @@
namespace scalfmm::interpolation
{
struct empty
{
};
template<typename ValueType, std::size_t Dimension, typename FarFieldMatrixKernel>
struct chebyshev_interpolator
: public interpolator<chebyshev_interpolator<ValueType, Dimension, FarFieldMatrixKernel>>
......@@ -68,58 +71,91 @@ namespace scalfmm::interpolation
using k_tensor_type = xt::xarray<value_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_inner_type = empty;
// using buffer_shape_type = xt::xshape<kn>;
// using buffer_type = xt::xtensor_fixed<buffer_inner_type, buffer_shape_type>;
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 operators = (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));
}
/**
* Chebyshev polynomials of first kind \f$ T_n(x) = \cos(n \arccos(x)) \f$
* 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 /*order*/, std::size_t n) const
[[nodiscard]] inline auto polynomials_impl(ComputationType x, std::size_t n) const
{
using value_type = ComputationType;
// assert(xsimd::any(xsimd::abs(x) - 1. < 10. * std::numeric_limits<value_type>::epsilon()));
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 value_type(cos(n * acos(x)));
return coeff * L;
}
/**
* 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] x coordinate in [-1,1]
* @param[in] order not used (just for comptibility)
* 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 /*order*/, std::size_t n) const
[[nodiscard]] inline auto derivative_impl(ComputationType x, std::size_t n) const
{
using value_type = ComputationType;
// assert(xsimd::any(xsimd::abs(x) - 1. < 10. * std::numeric_limits<value_type>::epsilon()));
simd::compare(x);
return value_type(n * (sin(n * acos(x))) / sqrt(1 - x * 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 {}
......@@ -134,16 +170,7 @@ namespace scalfmm::interpolation
return xt::fftw::rfftn<dimension>(C);
}
[[nodiscard]] inline auto buffer_initialization_impl() const -> buffer_type
{
struct empty
{
};
return empty{};
// 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.));
}
[[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
......@@ -208,15 +235,68 @@ namespace scalfmm::interpolation
// 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)));
// 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>
......@@ -226,12 +306,12 @@ namespace scalfmm::interpolation
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<svalue_type>;
// using buffer_inner_type = xt::xarray<std::complex<value_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_type = xt::xtensor_fixed<buffer_inner_type, xt::xshape<far_field_matrix_kernel_type::kn>>;
using buffer_inner_type = empty;
using buffer_type = empty;
};
} // namespace scalfmm::interpolation
......
......@@ -96,6 +96,7 @@ namespace scalfmm::interpolation
, m_nnodes(meta::pow(order, dimension))
, m_m2l_interactions(math::pow(max_number_of_cell, dimension))
{
this->derived_cast().set_private_members_impl();
set_child_parent_interpolator();
generate_interactions_matrixes(
homogeneity_tag == matrix_kernels::homogeneity::non_homogenous ? root_cell_width : value_type(2.),
......@@ -114,6 +115,7 @@ namespace scalfmm::interpolation
, m_nnodes(meta::pow(order, dimension))
, m_m2l_interactions(math::pow(max_number_of_cell, dimension))
{
this->derived_cast().set_private_members_impl();
set_child_parent_interpolator();
generate_interactions_matrixes(
homogeneity_tag == matrix_kernels::homogeneity::non_homogenous ? root_cell_width : value_type(2.),
......@@ -154,15 +156,15 @@ namespace scalfmm::interpolation
[[nodiscard]] inline auto roots(std::size_t order) const { return this->derived_cast().roots_impl(order); }
template<typename ComputationType>
[[nodiscard]] inline auto polynomials(ComputationType x, std::size_t order, std::size_t n) const
[[nodiscard]] inline auto polynomials(ComputationType x, std::size_t /*order*/, std::size_t n) const
{
return this->derived_cast().template polynomials_impl<ComputationType>(x, order, n);
return this->derived_cast().template polynomials_impl<ComputationType>(x, n);
}
template<typename ComputationType>
[[nodiscard]] inline auto derivative(ComputationType x, std::size_t order, std::size_t n) const
{
return this->derived_cast().template derivative_impl<ComputationType>(x, order, n);
return this->derived_cast().template derivative_impl<ComputationType>(x, n);
}
[[nodiscard]] inline auto separation_criterion() const noexcept -> std::size_t
......@@ -200,7 +202,7 @@ namespace scalfmm::interpolation
{
// rool level is 0 (i.e the simulation box) first level of cell is level 1
// we start the indexing of matrixes at 0, ence the -1 on the cell level
level = tree_level-2;
level = tree_level - 2;
}
return this->derived_cast().template apply_m2l_impl<Cell>(
source_cell, target_cell, m_interactions_matrixes.at(level * m_m2l_interactions + neighbor_idx),
......
......@@ -11,6 +11,7 @@
#include "scalfmm/interpolation/interpolator.hpp"
#include "scalfmm/interpolation/mapping.hpp"
#include "scalfmm/interpolation/traits.hpp"
#include "scalfmm/matrix_kernels/mk_common.hpp"
#include "scalfmm/meta/traits.hpp"
#include "scalfmm/meta/utils.hpp"
#include "scalfmm/simd/memory.hpp"
......@@ -18,7 +19,6 @@
#include "scalfmm/tools/colorized.hpp"
#include "scalfmm/utils/math.hpp"
#include "scalfmm/utils/tensor.hpp"
#include "scalfmm/matrix_kernels/mk_common.hpp"
#include <algorithm>
#include <array>
......@@ -51,7 +51,6 @@
#include <xtl/xclosure.hpp>
#include <xtl/xcomplex.hpp>
namespace scalfmm::interpolation
{
// uniform_interpolator is a CRTP based class
......@@ -88,7 +87,7 @@ namespace scalfmm::interpolation
}
template<typename ComputationType>
[[nodiscard]] inline auto polynomials_impl(ComputationType x, std::size_t order, std::size_t n) const
[[nodiscard]] inline auto polynomials_impl(ComputationType x, std::size_t n) const
{
using value_type = ComputationType;
// assert(xsimd::any(xsimd::abs(x) - 1. < 10. * std::numeric_limits<value_type>::epsilon()));
......@@ -96,7 +95,7 @@ namespace scalfmm::interpolation
const value_type two_const{2.0};
simd::compare(x);
auto order = this->order();
// Specific precomputation of scale factor
// in order to minimize round-off errors
// NB: scale factor could be hardcoded (just as the roots)
......@@ -137,12 +136,13 @@ namespace scalfmm::interpolation
}
template<typename ComputationType>
[[nodiscard]] inline auto derivative_impl(ComputationType x, std::size_t order, std::size_t n) const
[[nodiscard]] inline auto derivative_impl(ComputationType x, 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};
auto order{this->order()};
simd::compare(x);
// optimized variant
......@@ -245,7 +245,7 @@ namespace scalfmm::interpolation
// meta loop on kn
meta::repeat_indexed(
[&k, &products, &transformed_multipoles, m](std::size_t n, auto const& s_f) {
tensor::product(products.at(n), transformed_multipoles.at(m), k.at(n,m), s_f);
tensor::product(products.at(n), transformed_multipoles.at(m), k.at(n, m), s_f);
},
scale_factor);
}
......@@ -271,14 +271,14 @@ namespace scalfmm::interpolation
target_expansion);
}
//This function generates the circulent tensors for generating interaction matrixes.
// This function generates the circulent tensors for generating interaction matrixes.
template<typename TensorViewX, typename TensorViewY>
[[nodiscard]] inline auto generate_matrix_k(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
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{};
//generate the circulent tensor
// generate the circulent tensor
auto C = build_c(std::forward<TensorViewX>(X), std::forward<TensorViewY>(Y), order, far_field, n, m);
// return the fft
return xt::fftw::rfftn<dimension>(C);
......@@ -296,15 +296,15 @@ namespace scalfmm::interpolation
const value_type quarter{0.25};
// we get the half width to scale the roots
if constexpr (base_type::homogeneity_tag == matrix_kernels::homogeneity::non_homogenous)
if constexpr(base_type::homogeneity_tag == matrix_kernels::homogeneity::non_homogenous)
{
// (tree_heigh = 4 -> [0-3]) the bottom cells and leaves have the same symbolic level !
// but we remove level 0 (the root ie simulaton box) and level 1.
number_of_level = tree_height-2;
current_width = width*quarter;
number_of_level = tree_height - 2;
current_width = width * quarter;
}
value_type half_width{current_width*half};
value_type half_width{current_width * half};
// we need to keep the tensorial view
// let the same view goes down to tensorial of point
// X and Y are Td tensors storing point.
......@@ -317,11 +317,12 @@ namespace scalfmm::interpolation
// resizing and initializing the vector
// homogenous -> only one level
// non_homogenous -> we generate interaction matrices for each tree level processed by the m2l operator
interactions_matrixes.resize(number_of_interactions*number_of_level, interaction_matrix_type(typename interaction_matrix_type::shape_type{},
k_tensor_type(std::vector(dimension, this->order())),
xt::layout_type::row_major));
interactions_matrixes.resize(number_of_interactions * number_of_level,
interaction_matrix_type(typename interaction_matrix_type::shape_type{},
k_tensor_type(std::vector(dimension, this->order())),
xt::layout_type::row_major));
// loop on levels
for(std::size_t l{0}; l<number_of_level; ++l)
for(std::size_t l{0}; l < number_of_level; ++l)
{
// homogenous -> X is the [-1,1] box reference
// non_homogenous -> X is the [-cell_width_at_level/2, cell_width_at_level/2]
......@@ -335,11 +336,11 @@ namespace scalfmm::interpolation
// get an array of points
auto X_points = xt::xarray<container::point<value_type, dimension>>::from_shape(std::get<0>(X).shape());
//we flatten the grids
// we flatten the grids
auto X_flatten_views = std::apply(
[](auto&&... xs) { return std::make_tuple(xt::flatten(std::forward<decltype(xs)>(xs))...); }, X);
//we flatten the tensor of points
// we flatten the tensor of points
auto X_points_flatten_views = xt::flatten(X_points);
// here, we reconstruct the points for the grids.
......@@ -355,7 +356,8 @@ namespace scalfmm::interpolation
// lambda for generation the matrixes corresponding to one interaction
std::size_t flat_idx{0};
auto generate_all_interactions = [this, &interactions_matrixes, &flat_idx, &X_points, current_width, number_of_interactions, l](auto... is) {
auto generate_all_interactions = [this, &interactions_matrixes, &flat_idx, &X_points, current_width,
number_of_interactions, l](auto... is) {
if(((std::abs(is) > this->separation_criterion()) || ...))
{
// constructing centers to generate Y
......@@ -370,7 +372,7 @@ namespace scalfmm::interpolation
Y_points = X_points + centers;
// we get the ref of interaction matrices to generate
auto& nm_fc_tensor = interactions_matrixes.at(l*number_of_interactions + flat_idx);
auto& nm_fc_tensor = interactions_matrixes.at(l * number_of_interactions + flat_idx);
// and we generate each matrix needed for the product (kn*km matrices)
for(std::size_t n = 0; n < kn; ++n)
{
......@@ -384,7 +386,7 @@ namespace scalfmm::interpolation
++flat_idx;
};
//loop range [-3,4[
// loop range [-3,4[
std::array<int, dimension> starts{};
std::array<int, dimension> stops{};
starts.fill(-3);
......@@ -398,7 +400,7 @@ namespace scalfmm::interpolation
half_width *= half;
}
}
inline void set_private_members_impl() {}
};
// traits class to register types inside the interpolator generic class
......
......@@ -6,6 +6,7 @@ set(source_tests_files
container/particle.cpp
tree/interaction_list.cpp
operators/l2p.cpp
operators/interp.cpp
fmm/test_dimension.cpp
fmm/count_kernel.cpp
fmm/test_non_homogenous.cpp
......
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