Commit 3c0ef569 authored by ESTERIE Pierre's avatar ESTERIE Pierre

Add potentials and multipoles calculation in p2m operator

parent 86b0d024
......@@ -6,6 +6,7 @@
#include <iostream>
#include <utility>
#include <random>
#include <tuple>
#include <scalfmm/container/particle.hpp>
#include <scalfmm/container/particle_container.hpp>
......@@ -66,14 +67,47 @@ int main()
// Construct leaf
leaf_type l(std::move(particles_source), center, width);
multipole_type multipoles{scalfmm::meta::pow(order,particle_dim)};
start = std::chrono::high_resolution_clock::now();
scalfmm::operators::apply_p2m(s, order, l, multipoles);
end = std::chrono::high_resolution_clock::now();
multipole_type multipoles(scalfmm::meta::pow(order,particle_dim), std::make_tuple(0.,0.,0.));
std::vector<std::chrono::microseconds> t;
for(std::size_t exp = 0; exp < 50; ++exp)
{
//multipoles.clear();
start = std::chrono::high_resolution_clock::now();
scalfmm::operators::apply_p2m(s, order, l, multipoles);
end = std::chrono::high_resolution_clock::now();
t.push_back(std::chrono::duration_cast<std::chrono::microseconds>(end - start));
if(exp == 0)
{
std::ofstream file("dump.txt");
for(int i = 0; i < multipoles.size(); ++i)
{
auto it = std::begin(multipoles);
file << std::fixed << std::setprecision(13) << std::get<0>(*it) << ' ';
it++;
}
for(int i = 0; i < multipoles.size(); ++i)
{
auto it = std::begin(multipoles);
file << std::fixed << std::setprecision(13) << std::get<1>(*it) << ' ';
it++;
}
for(int i = 0; i < multipoles.size(); ++i)
{
auto it = std::begin(multipoles);
file << std::fixed << std::setprecision(13) << std::get<2>(*it) << ' ';
it++;
}
file.close();
}
}
std::sort(t.begin(),t.end());
std::cout << colors::yellow << "p2m took... : "
<< std::chrono::duration_cast<std::chrono::microseconds>(end - start).count() << "us\n"
<< t[t.size()/2].count() << "us\n"
<< colors::reset;
return 0;
......
......@@ -49,6 +49,12 @@ namespace scalfmm::simd
return std::make_tuple(xsimd::load_aligned(std::get<Particle::dimension+Is>(particle_ptrs))...);
}
template<typename TupleAttributes, typename TuplePtr, std::size_t... Is>
constexpr inline void store_attributes(TuplePtr const& attributes_ptrs, TupleAttributes const& src,
std::index_sequence<Is...> s)
{
noop_f{(xsimd::store_aligned(std::get<Is>(attributes_ptrs), std::get<Is>(src)), 0)...};
}
template<typename F, typename Container, typename... Types, std::size_t... Is>
constexpr inline auto apply_f(std::index_sequence<Is...> s, F&& f, Container&& input, Types&&... ts)
......@@ -77,16 +83,16 @@ namespace scalfmm::simd
return impl::load_attributes<Particle>(particle_ptrs, std::make_index_sequence<Particle::attributes_size>{});
}
template<typename Particle, typename Position, typename TuplePtr>
template<typename Position, typename TuplePtr>
constexpr inline void store_position(TuplePtr& particle_ptrs, Position const& src)
{
impl::store_position(particle_ptrs, src, std::make_index_sequence<Particle::dimension>{});
impl::store_position(particle_ptrs, src, std::make_index_sequence<Position::dimension>{});
}
template<typename Particle, std::size_t Index, typename TuplePtr, typename Attribute>
constexpr inline void store_attribute(TuplePtr& particle_ptrs, Attribute const& src)
template<typename TupleAttributes, typename TuplePtr>
constexpr inline void store_attributes(TuplePtr& attributes_ptrs, TupleAttributes const& src)
{
xsimd::store_aligned(std::get<Particle::dimension>(particle_ptrs), src);
impl::store_attributes(attributes_ptrs, src, std::make_index_sequence<std::tuple_size<TupleAttributes>::value>{});
}
template<std::size_t Size, typename F, typename Container, typename... Types>
......@@ -103,6 +109,25 @@ namespace scalfmm::simd
}
} // namespace scalfmm::simd
namespace scalfmm::utils
{
namespace impl
{
template<typename T, typename U, std::size_t... Is>
inline constexpr auto tuple_update(T&& lhs, U&& rhs, std::index_sequence<Is...> s)
{
simd::impl::noop_f{((std::get<Is>(std::forward<T>(lhs))+=std::get<Is>(std::forward<U>(rhs))),0)...};
}
}
template<typename T, typename U>
inline constexpr void tuple_update(T&& lhs, U&& rhs)
{
return impl::tuple_update(std::forward<T>(lhs), std::forward<U>(rhs), std::make_index_sequence<std::tuple_size<T>::value>{});
}
}
namespace scalfmm::operators
{
template<typename Interpolator, typename Leaf, typename Multipoles>
......@@ -155,7 +180,7 @@ namespace scalfmm::operators
// Assembling S
std::size_t idx = 0;
auto construct_s = [/*&part,*/ &S, &poly_of_part, &idx, inc](auto&... current_indices) {
auto construct_s = [&S, &poly_of_part, &idx, inc](auto&... current_indices) {
//meta::print(std::cout, current_indices...);
//std::cout << '\n';
auto s_simd = simd::generate_s<position_type::dimension>(poly_of_part, {{current_indices...}});
......@@ -169,8 +194,39 @@ namespace scalfmm::operators
meta::looper<position_type::dimension>()(construct_s, stops);
auto simd_potentials = simd::load_attributes<particle_type>(particle_iterator);
//std::cout << std::get<0>(simd_potentials) << '\n';
auto update_potential = [](auto pot, const value_type* current_s)
{
auto s_to_apply = xsimd::load_aligned(current_s);
return s_to_apply*pot;
};
std::size_t vec_idx = 0;
for(std::size_t m = 0; m < expansion.size(); ++m)
{
auto s_to_apply = xsimd::load_aligned(&S[vec_idx]);
auto simd_potentials = simd::load_attributes<particle_type>(particle_iterator);
auto multipoles_iterator = expansion.begin();
auto update_potential = [&s_to_apply](auto... pot)
{
return std::make_tuple(s_to_apply*pot...);
};
simd_potentials = std::apply(update_potential, simd_potentials);
auto reduce_potentials = [](auto... pot)
{
return std::make_tuple(xsimd::hadd(pot)...);
};
auto weight_to_apply = std::apply(reduce_potentials, simd_potentials);
utils::tuple_update(*multipoles_iterator,weight_to_apply);
multipoles_iterator++;
vec_idx += inc;
}
poly_of_part.clear();
S.clear();
......
Markdown is supported
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