Attention une mise à jour du service Gitlab va être effectuée le mardi 30 novembre entre 17h30 et 18h00. Cette mise à jour va générer une interruption du service dont nous ne maîtrisons pas complètement la durée mais qui ne devrait pas excéder quelques minutes. Cette mise à jour intermédiaire en version 14.0.12 nous permettra de rapidement pouvoir mettre à votre disposition une version plus récente.

Commit 8f71e826 authored by ESTERIE Pierre's avatar ESTERIE Pierre
Browse files

some new tree construction + time loop example + exportVTK from tree

parent d1f0e72b
......@@ -22,6 +22,7 @@ set(source_tests_files
test_l2p.cpp
test-build-let.cpp
trace-tree.cpp
test-time-loop.cpp
)
#if(SCALFMM_USE_MPI)
# set(source_tests_files ${source_tests_files}
......

#include "parameters.hpp"
#include "scalfmm/algorithms/common.hpp"
#include "scalfmm/container/particle.hpp"
#include "scalfmm/container/particle_container.hpp"
#include "scalfmm/container/point.hpp"
#include "scalfmm/interpolation/uniform.hpp"
#include "scalfmm/operators/fmmOperators.hpp"
#include "scalfmm/matrix_kernels/laplace.hpp"
#include "scalfmm/matrix_kernels/scalar_kernels.hpp"
//#include "scalfmm/operators/generic/l2p.hpp"
#include "scalfmm/algorithms/sequential/sequential.hpp"
#include "scalfmm/tools/colorized.hpp"
#include "scalfmm/tools/tikz_writer.hpp"
#include "scalfmm/tree/box.hpp"
#include "scalfmm/tree/cell.hpp"
#include "scalfmm/tree/group_tree.hpp"
#include "scalfmm/tree/leaf.hpp"
#include "scalfmm/utils/sort.hpp"
#include "scalfmm/tools/fma_loader.hpp"
#include "scalfmm/tools/vtk_writer.hpp"
#include "scalfmm/utils/accurater.hpp"
#include "scalfmm/utils/timer.hpp"
#include <algorithm>
#include <array>
#include <bits/c++config.h>
#include <chrono>
#include <cstdio>
#include <inria/tcli/tcli.hpp>
#include <iostream>
#include <sstream>
#include <string>
#include <sys/types.h>
#include <thread>
#include <tuple>
#include <unistd.h>
#include <utility>
#include <vector>
// example:
// ./examples/RelWithDebInfo/test-time-loop --tree-height 4 -gs 2 --order 4 --dimension 2...
namespace local_args
{
struct tree_height : inria::tcli::required_tag
{
inria::tcli::str_vec flags = {"--tree-height", "-th"};
const char* description = "Tree height (or initial height in case of an adaptive tree).";
using type = std::size_t;
type def = 2;
};
struct order : inria::tcli::required_tag
{
inria::tcli::str_vec flags = {"--order", "-o"};
const char* description = "Precision setting.";
using type = std::size_t;
type def = 3;
};
struct thread_count
{
inria::tcli::str_vec flags = {"--threads", "-t"};
const char* description = "Maximum thread count to be used.";
using type = std::size_t;
type def = 1;
};
struct input_file : inria::tcli::required_tag
{
inria::tcli::str_vec flags = {"--input-file", "-fin"};
const char* description = "Input filename (.fma or .bfma).";
using type = std::string;
};
struct output_file
{
inria::tcli::str_vec flags = {"--output-file", "-fout"};
const char* description = "Output particle file (with extension .fma (ascii) or bfma (binary).";
using type = std::string;
};
struct visu_file
{
inria::tcli::str_vec flags = {"--visu-file", "-vf"};
const char* description = "Output VTK file.";
using type = std::string;
};
struct block_size
{
inria::tcli::str_vec flags = {"--group-size", "-gs"};
const char* description = "Group tree chunk size.";
using type = int;
type def = 250;
};
struct dimension
{
inria::tcli::str_vec flags = {"--dimension", "--d"};
const char* description = "Dimension : \n 2 for dimension 2, 3 for dimension 3";
using type = std::size_t;
type def = 1;
};
struct time_steps
{
inria::tcli::str_vec flags = {"--time-step", "--ts"};
const char* description = "Number of time steps";
using type = std::size_t;
type def = 1;
};
} // namespace local_args
template<std::size_t Dimension, typename ContainerType, typename ValueType>
void read_data(const std::string& filename, ContainerType*& container,
scalfmm::container::point<ValueType, Dimension>& Centre, ValueType& width)
{
using particle_type = typename ContainerType::value_type;
scalfmm::tools::FFmaGenericLoader<ValueType, Dimension> loader(filename);
const int number_of_particles = loader.getNumberOfParticles();
width = loader.getBoxWidth();
Centre = loader.getBoxCenter();
auto nb_val_to_red_per_part = loader.getNbRecordPerline();
double* values_to_read = new double[nb_val_to_red_per_part]{0};
container = new ContainerType(number_of_particles);
std::cout << "number_of_particles " << number_of_particles << std::endl;
for(std::size_t idx = 0; idx < number_of_particles; ++idx)
{
loader.fillParticle(values_to_read, nb_val_to_red_per_part);
particle_type p;
std::size_t ii{0};
for(auto& e: p.position())
{
e = values_to_read[ii++];
}
for(auto& e: p.inputs())
{
e = ValueType(1.0);
}
// p.variables(values_to_read[ii++], idx);
container->insert_particle(idx, p);
}
}
template<std::size_t Dimension, typename FmmOperatorType, typename... Parameters>
auto run(inria::tcli::parser<Parameters...> const& parser) -> int
{
std::cout << scalfmm::colors::blue << "Entering tree test...\n" << scalfmm::colors::reset;
// timer
scalfmm::utils::timer time{};
using value_type = double;
// ---------------------------------------
using near_matrix_kernel_type = typename FmmOperatorType::near_field_type::matrix_kernel_type;
using far_field_type = typename FmmOperatorType::far_field_type;
using interpolator_type = typename far_field_type::approximation_type;
using far_matrix_kernel_type = typename interpolator_type::matrix_kernel_type;
// ---------------------------------------
// The matrix kernel
static constexpr std::size_t nb_inputs_near{near_matrix_kernel_type::km};
static constexpr std::size_t nb_outputs_near{near_matrix_kernel_type::kn};
static constexpr std::size_t nb_inputs_far{far_matrix_kernel_type::km};
static constexpr std::size_t nb_outputs_far{far_matrix_kernel_type::kn};
// ---------------------------------------
using particle_type = scalfmm::container::particle<value_type, Dimension, value_type, nb_inputs_near, value_type,
nb_outputs_near, std::size_t>;
using container_type = scalfmm::container::particle_container<particle_type>;
using position_type = typename particle_type::position_type;
using cell_type = scalfmm::component::cell<typename interpolator_type::storage_type>;
using leaf_type = scalfmm::component::leaf<particle_type>;
using box_type = scalfmm::component::box<position_type>;
using group_tree_type = scalfmm::component::group_tree<cell_type, leaf_type, box_type>;
// ---------------------------------------
// parameters
const auto tree_height{parser.template get<local_args::tree_height>()};
const auto group_size{parser.template get<local_args::block_size>()};
const auto order{parser.template get<local_args::order>()};
const auto n_time_steps = parser.template get<local_args::time_steps>();
const std::string input_file(parser.template get<local_args::input_file>());
const std::string output_file(parser.template get<local_args::output_file>());
const std::string visu_file(parser.template get<local_args::visu_file>());
std::cout << scalfmm::colors::blue << "<params> Tree height : " << tree_height << scalfmm::colors::reset << '\n';
std::cout << scalfmm::colors::blue << "<params> Group Size : " << group_size << scalfmm::colors::reset << '\n';
std::cout << scalfmm::colors::blue << "<params> Runtime order : " << order << scalfmm::colors::reset << '\n';
std::cout << scalfmm::colors::blue << "<params> Time Steps : " << n_time_steps << scalfmm::colors::reset << '\n';
// ---------------------------------------
// Open particle file
std::cout << scalfmm::colors::green << "Creating & Inserting particles ...\n" << scalfmm::colors::reset;
scalfmm::container::point<value_type, Dimension> box_center{};
value_type box_width{};
container_type* container{};
if(input_file.empty())
{
std::cerr << scalfmm::colors::red << "input file is empty !\n";
return 0;
}
else
{
time.tic();
read_data(input_file, container, box_center, box_width);
time.tac();
}
std::cout << scalfmm::colors::green << "... Done.\n" << scalfmm::colors::reset;
std::cout << scalfmm::colors::green << "Box center = " << box_center << " box width = " << box_width
<< scalfmm::colors::reset << '\n';
std::cout << scalfmm::colors::yellow << "Container loaded in " << time.elapsed() << "ms\n"
<< scalfmm::colors::reset;
// ---------------------------------------
// creating group_tree
time.tic();
box_type box(box_width, box_center);
group_tree_type tree(tree_height, order, box, group_size, group_size, *container);
time.tac();
std::cout << scalfmm::colors::yellow << "Group tree created in " << time.elapsed() << "ms\n"
<< scalfmm::colors::reset;
// ---------------------------------------
// creating interpolator
time.tic();
far_matrix_kernel_type mk_far{};
interpolator_type interpolator(mk_far, order, static_cast<std::size_t>(tree_height), box.width(0));
near_matrix_kernel_type mk_near{};
typename FmmOperatorType::near_field_type near_field(mk_near);
typename FmmOperatorType::far_field_type far_field(interpolator);
FmmOperatorType fmm_operator(near_field, far_field);
time.tac();
std::cout << scalfmm::colors::yellow << "Kernel and Interp created in " << time.elapsed() << "ms\n"
<< scalfmm::colors::reset;
// ---------------------------------------
// start algorithm for step 0
auto operator_to_proceed = scalfmm::algorithms::all;
scalfmm::algorithms::sequential(tree, std::move(fmm_operator), operator_to_proceed);
// ---------------------------------------
// time loop
for(std::size_t steps{1}; steps < n_time_steps; ++steps)
{
// ---
// compute the forces
// ---
// update particles
// ---
// compute the corners and set new box
// ---
// compute new morton indices
// ---
// move particles according to new morton indices
// ---
// redistribute particle in leaves
// ---
// build tree levels
// ---
}
// ---------------------------------------
// write output fma
if(!output_file.empty())
{
std::cout << "Write outputs in " << output_file << std::endl;
scalfmm::tools::FFmaGenericWriter<double> writer(output_file);
writer.writeDataFromTree(tree, container->size());
}
// ---------------------------------------
// write output fma
if(!visu_file.empty())
{
std::cout << "Write VTK file in " << visu_file << std::endl;
scalfmm::tools::io::exportVTKxml(visu_file, tree, container->size());
}
delete container;
return 0;
}
auto main([[maybe_unused]] int argc, [[maybe_unused]] char* argv[]) -> int
{
//
// Parameter handling
auto parser =
inria::tcli::make_parser(inria::tcli::help{}, local_args::tree_height{}, local_args::order{},
local_args::input_file{}, local_args::output_file{}, local_args::block_size{},
local_args::dimension{}, local_args::time_steps{}, local_args::visu_file{});
parser.parse(argc, argv);
// Getting command line parameters
const std::size_t dimension = parser.get<local_args::dimension>();
switch(dimension)
{
case 2:
{
constexpr std::size_t dim = 2;
using matrix_kernel_type = scalfmm::matrix_kernels::others::one_over_r2;
using near_field_type = scalfmm::operators::near_field_operator<matrix_kernel_type>;
//
using interpolator_type = scalfmm::interpolation::uniform_interpolator<double, dim, matrix_kernel_type>;
using far_field_type = scalfmm::operators::far_field_operator<interpolator_type>;
run<dim, scalfmm::operators::fmm_operators<near_field_type, far_field_type>>(parser);
break;
}
case 3:
{
constexpr std::size_t dim = 3;
using matrix_kernel_type = scalfmm::matrix_kernels::others::one_over_r2;
// using matrix_kernel_type = scalfmm::matrix_kernels::laplace::one_over_r;
using near_field_type = scalfmm::operators::near_field_operator<matrix_kernel_type>;
//
using interpolator_type = scalfmm::interpolation::uniform_interpolator<double, dim, matrix_kernel_type>;
using far_field_type = scalfmm::operators::far_field_operator<interpolator_type>;
run<dim, scalfmm::operators::fmm_operators<near_field_type, far_field_type>>(parser);
break;
}
default:
std::cout << "check 1/r^2 Kernel for dimension 2 and 3. Value is \n"
<< " 0 for dimension 2 "
<< " 1 for dimension 3 " << std::endl;
break;
}
}
......@@ -5,9 +5,12 @@
#ifndef SCALFMM_TOOLS_VTK_WRITER_HPP
#define SCALFMM_TOOLS_VTK_WRITER_HPP
#include <bits/c++config.h>
#include <fstream>
#include <iostream>
#include <string>
#include "scalfmm/meta/utils.hpp"
#include "scalfmm/tree/for_each.hpp"
namespace scalfmm::tools::io
{
......@@ -26,6 +29,136 @@ namespace scalfmm::tools::io
return header;
}
template<class TreeType>
void exportVTKxml(std::string const& filename, TreeType const& tree, std::size_t npart)
{
using value_type = typename TreeType::leaf_type::value_type;
using particle_type = typename TreeType::particle_type;
std::vector<particle_type> particles(npart);
std::size_t pos{0};
scalfmm::component::for_each_leaf(std::cbegin(tree), std::cend(tree), [&pos,&particles](auto& leaf) {
const auto container = leaf.cparticles();
for(std::size_t idx = 0; idx < leaf.cparticles().size(); ++idx)
{
particles[pos++] = leaf.cparticles().particle(idx);
}
});
if(TreeType::dimension > 3)
{
std::cerr << " exportVTKxml works only for dimension <=3 and the "
"results are in space 3"
<< std::endl;
return;
}
std::ofstream VTKfile(filename);
std::size_t j = 0;
std::size_t stride = meta::tuple_size_v<particle_type>;
int nb_output_values = particle_type::outputs_size;
int nb_input_values = particle_type::inputs_size;
std::cout << "Write vtk format in " << filename << std::endl;
std::cout << " dim=" << TreeType::dimension << " stride =" << stride << " nb_input_values " << nb_input_values
<< " nb_output_values " << nb_output_values << " N " << npart << std::endl;
VTKfile << "<?xml version=\"1.0\"?>" << std::endl
<< "<VTKFile type=\"PolyData\" version=\"0.1\" "
"byte_order=\"LittleEndian\"> "
<< std::endl
<< "<PolyData>" << std::endl
<< "<Piece NumberOfPoints=\" " << particles.size() << " \" NumberOfVerts=\" " << particles.size()
<< " \" NumberOfLines=\" 0\" NumberOfStrips=\"0\" NumberOfPolys=\"0\">" << std::endl
<< "<Points>" << std::endl
<< "<DataArray type=\"Float64\" NumberOfComponents=\"3" /*<< dimension*/
<< "\" "
"format=\"ascii\"> "
<< std::endl;
j = 0;
for(std::size_t i = 0; i < particles.size(); ++i)
{
for(int k = 0; k < TreeType::dimension; ++k)
{
VTKfile << particles[i].position(k) << " ";
}
for(int i = 0; i < 3 - TreeType::dimension; ++i)
{
VTKfile << "0.0 ";
}
}
if(nb_input_values + nb_output_values > 0)
{
std::string header("Scalars=\"");
if(nb_input_values > 0)
{
for(int i = 0; i < nb_input_values; ++i)
{
header += "intputValue" + std::to_string(i) + " ";
}
for(int i = 0; i < nb_output_values; ++i)
{
header += "outputValue" + std::to_string(i) + " ";
}
header += "\" ";
}
VTKfile << std::endl
<< "</DataArray> " << std::endl
<< "</Points> " << std::endl
<< "<PointData " << header << " > " << std::endl;
for(int k = 0; k < nb_input_values; ++k)
{
VTKfile << "<DataArray type=\"Float64\" Name=\"inputValue" + std::to_string(k) +
"\" "
"format=\"ascii\">"
<< std::endl;
j = 0;
for(std::size_t i = 0; i < particles.size(); ++i)
{
VTKfile << particles[i].inputs(k) << " ";
}
VTKfile << std::endl << "</DataArray>" << std::endl;
}
for(int k = 0; k < nb_output_values; ++k)
{
VTKfile << "<DataArray type=\"Float64\" Name=\"outputValue" + std::to_string(k - nb_input_values) +
"\" "
"format=\"ascii\">"
<< std::endl;
j = 0;
for(std::size_t i = 0; i < particles.size(); ++i)
{
VTKfile << particles[i].outputs(k) << " ";
}
VTKfile << std::endl << "</DataArray>" << std::endl;
}
VTKfile << " </PointData>" << std::endl
<< " <CellData>"
<< " </CellData>" << std::endl;
}
VTKfile << " <Verts>" << std::endl
<< " <DataArray type=\"Int32\" Name=\"connectivity\" "
"format=\"ascii\">"
<< std::endl;
for(std::size_t i = 0; i < particles.size(); ++i)
{
VTKfile << i << " ";
}
VTKfile << std::endl
<< "</DataArray>" << std::endl
<< "<DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << std::endl;
for(std::size_t i = 1; i < particles.size() + 1; ++i)
{
VTKfile << i << " ";
}
VTKfile << std::endl
<< "</DataArray>" << std::endl
<< " </Verts>" << std::endl
<< "<Lines></Lines>" << std::endl
<< "<Strips></Strips>" << std::endl
<< "<Polys></Polys>" << std::endl
<< "</Piece>" << std::endl
<< "</PolyData>" << std::endl
<< "</VTKFile>" << std::endl;
}
//! \fn void exportVTKxml(std::ofstream& file, const FReal particles, const
//! std::size_t N )
......
......@@ -441,55 +441,6 @@ namespace scalfmm::component
}
}
/// @brief The function fills the particle in each leaf.
/// Its uses the source index in the tuple of indices to get the
/// particle from the source container
///
/// @tparam VectorOfTuples : a vector of tuples of indices of size number_of_leaves
/// @tparam ParticleContainer : the particle container
/// @param tuple_of_indexes : a vector holding indices for permutation of size number_of_leaves
/// @param particle_container : the container storing all the particles.
///
/// @return
template<typename VectorOfTuples, typename ParticleContainer>
auto fill_leaves_with_particles(VectorOfTuples const& tuple_of_indexes,
ParticleContainer const& particle_container) -> void
{
using proxy_type = typename particle_type::proxy_type;
using outputs_value_type = typename particle_type::outputs_value_type;
auto begin_container = std::begin(particle_container);
std::size_t leaf_index{0};
// loop on groups
for(auto pg: m_group_of_leaf)
{
// loop on leaves
for(auto& leaf: pg->block())
{
// get the leaf container
auto leaf_container_begin = std::begin(leaf.particles());
// copy the particle in the leaf
for(std::size_t index_part = 0; index_part < leaf.size(); ++index_part)
{
// get the source index in the source container
auto source_index = std::get<1>(tuple_of_indexes.at(leaf_index++));
// jump to the index in the source container
auto jump_to_particle = begin_container;
std::advance(jump_to_particle, source_index);
// copy the particle
*leaf_container_begin = meta::as_tuple(*jump_to_particle);
proxy_type particle_in_leaf(*leaf_container_begin);
// set the outputs to zero
for(std::size_t ii{0}; ii < particle_type::outputs_size; ++ii)
{
particle_in_leaf.outputs(ii) = outputs_value_type(0.);
}
++leaf_container_begin;
}
}
}
}
public:
/// @brief Constructor of the group tree.
/// It initialized all levels with leaves and cells from the particle container passed.
......@@ -520,9 +471,6 @@ namespace scalfmm::component
"group_tree : Particles contain in leafs are not the same as the ones in the container passed "
"as argument to this constructor.");
std::cout << "Component per leaf group = " << m_number_of_leaves_per_group << '\n';
std::cout << "Component per cell group = " << m_number_of_cells_per_group << '\n';
auto number_of_particles{particle_container.size()};
const std::size_t leaf_level{m_tree_height - 1};
......@@ -539,32 +487,27 @@ namespace scalfmm::component
std::vector<std::size_t> vector_of_mortons(tuple_of_indexes.size());
std::transform(std::begin(tuple_of_indexes), std::end(tuple_of_indexes), std::begin(vector_of_mortons),
[](auto const& t) { return std::get<0>(t); });
// vector to store the number of particles per leaves
auto number_of_particles_per_leaves{get_leaves_distribution(vector_of_mortons)};
// we build the first level of group for the leaves
build_groups_of_leaves(vector_of_mortons);
// we construct the leaves in each group
build_leaves_in_groups(vector_of_mortons, number_of_particles_per_leaves, box);
// construct all levels
construct(vector_of_mortons);
// then, we fill each leaf with its particle
fill_leaves_with_particles(tuple_of_indexes, particle_container);
// construct group of cells at leaf level
build_groups_of_cells_at_level(vector_of_mortons, leaf_level);
// construct cells in group of leaf level
build_cells_in_groups_at_level(vector_of_mortons, box, leaf_level);
}
// loop on levels
for(std::size_t level{m_tree_height - 2}; level >= m_top_level; --level)
{
// update vector_of_mortons for upper level
index::get_parent_morton_indices(vector_of_mortons, 0, dimension);
// construct group of cells at current level
build_groups_of_cells_at_level(vector_of_mortons, level);
// construct cells in group of current level
build_cells_in_groups_at_level(vector_of_mortons, box, level);
}
template<typename MortonType>
group_tree(std::size_t tree_height, std::size_t order, box_type const& box,
std::size_t number_of_leaves_per_group, std::size_t number_of_cells_per_group,
std::vector<MortonType> const& vector_of_mortons, int in_top_level = 2)