diff --git a/experimental/CMakeLists.txt b/experimental/CMakeLists.txt
index 1f0f9bdace1448f50e9b3756e4c3af2876003be5..64ac4c5bdf3f24d4bfcb21931cfbded1d707fec3 100644
--- a/experimental/CMakeLists.txt
+++ b/experimental/CMakeLists.txt
@@ -2,7 +2,7 @@
#
# Project Declaration
#--------------------
-project(scalfmm3 CXX)
+project(scalfmm3 CXX C)
# check if compiling into source directories
string(COMPARE EQUAL "${CMAKE_CURRENT_SOURCE_DIR}" "${CMAKE_CURRENT_BINARY_DIR}" insource)
if(insource)
diff --git a/experimental/cmake/dependencies.cmake b/experimental/cmake/dependencies.cmake
index 7b5e7d08a7358c99b010cb6aa8fe991b97f42970..c916bb5aee410684600ce5d78b599b5baf1e2648 100644
--- a/experimental/cmake/dependencies.cmake
+++ b/experimental/cmake/dependencies.cmake
@@ -8,3 +8,11 @@ include(cmake/dependencies/blas.cmake)
include(cmake/dependencies/fftw.cmake)
include(cmake/dependencies/starpu.cmake)
+if(SCALFMM_USE_CATALYST)
+ find_package(catalyst REQUIRED)
+ if(catalyst_FOUND)
+ message(STATUS "Catalyst found.")
+ target_link_libraries(${CMAKE_PROJECT_NAME} INTERFACE catalyst::catalyst catalyst::conduit_headers catalyst::conduit catalyst::blueprint)
+ endif()
+endif()
+
diff --git a/experimental/examples/catalyst_related.hpp b/experimental/examples/catalyst_related.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..67a79551d6412bc590f3b4c10a4927fb72eae897
--- /dev/null
+++ b/experimental/examples/catalyst_related.hpp
@@ -0,0 +1,89 @@
+// --------------------------------
+// See LICENCE file at project root
+// File : group_tree.hpp
+// --------------------------------
+#ifndef CATALYST_RELATED_HPP
+#define CATALYST_RELATED_HPP
+
+#include "scalfmm/meta/const_functions.hpp"
+#include "scalfmm/meta/utils.hpp"
+#include "scalfmm/utils/parallel_manager.hpp"
+#include <catalyst.h>
+#include <conduit.hpp>
+#include <conduit_blueprint_mesh.h>
+#include <conduit_blueprint_exports.h>
+#include <conduit_cpp_to_c.hpp>
+#include <inria/tcli/tcli.hpp>
+//#include <iostream>
+
+// catalyst
+namespace catalyst_adaptor
+{
+ template<typename... Parameters>
+ auto initialize(inria::tcli::parser<Parameters...> const& parser) -> void
+ {
+ parallel_manager m;
+ m.init();
+ // here goes scrpits of mpi init communicator
+ conduit::Node node;
+ //node["catalyst/scripts/script" + std::to_string(0)].set_string("/home/p13ro/dev/cpp/build/scalfmm/gcc-10.2/catalyst_pipeline.py");
+ node["catalyst/mpi_comm"].set(0);
+ catalyst_initialize(conduit::c_node(&node));
+ }
+
+ template<typename GroupTree>
+ auto execute(std::size_t step, GroupTree const& tree)
+ {
+ conduit::Node exec_params;
+
+ //// add time/cycle information
+ auto& state = exec_params["catalyst/state"];
+ state["timestep"].set(step);
+ state["time"].set(step);
+
+ // Add channels.
+ // We only have 1 channel here. Let's name it 'grid'.
+ //auto particles_channel = exec_params["catalyst/channels/particles"];
+ auto& box_channel = exec_params["catalyst/channels/box"];
+
+ // Since this example is using Conduit Mesh Blueprint to define the mesh,
+ // we set the channel's type to "mesh".
+ box_channel["type"].set("mesh");
+
+ const auto n_corners = scalfmm::meta::pow(2,GroupTree::dimension);
+
+ // now create the mesh.
+ auto& box_mesh = box_channel["data"];
+ box_mesh["coordsets/coords/type"].set("explicit");
+ box_mesh["coordsets/coords/values/x"].set(conduit::DataType::float64(n_corners));
+ box_mesh["coordsets/coords/values/y"].set(conduit::DataType::float64(n_corners));
+
+ auto const& box = tree.box();
+
+ conduit::float64* box_sim_x = box_mesh["coordsets/coords/values/x"].value();
+ conduit::float64* box_sim_y = box_mesh["coordsets/coords/values/x"].value();
+
+ for(std::size_t n{0}; n<n_corners; ++n)
+ {
+ box_sim_x[n] = scalfmm::meta::get<0>(box.corner(n));
+ box_sim_y[n] = scalfmm::meta::get<1>(box.corner(n));
+ std::cout << box.corner(n) << '\n';
+ }
+
+ // Next, add topology
+ box_mesh["topologies/mesh/type"].set("line");
+ box_mesh["topologies/mesh/coordset"].set("coords");
+ //box_mesh["topologies/mesh/elements/shape"].set("box");
+ box_mesh["topologies/mesh/elements/connectivity"].set_int32_vector({0,1,2,3});
+
+ catalyst_execute(conduit::c_node(&exec_params));
+ }
+
+ auto finalize() -> void
+ {
+ conduit::Node node;
+ catalyst_finalize(conduit::c_node(&node));
+ }
+}
+
+#endif // CATALYST_RELATED_HPP
diff --git a/experimental/examples/test-time-loop.cpp b/experimental/examples/test-time-loop.cpp
index 80f656f77d82da5d3eaacd2082a35e446ceea2b9..ebde79f67095a211c445fc3c057bf95bb9973f91 100644
--- a/experimental/examples/test-time-loop.cpp
+++ b/experimental/examples/test-time-loop.cpp
@@ -1,25 +1,20 @@
-
-#include "parameters.hpp"
-#include "scalfmm/algorithms/common.hpp"
-#include "scalfmm/container/particle.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/operators/laplace_operators.hpp"
#include "scalfmm/matrix_kernels/laplace.hpp"
#include "scalfmm/matrix_kernels/scalar_kernels.hpp"
-//#include "scalfmm/operators/generic/l2p.hpp"
+#include "scalfmm/algorithms/common.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/for_each.hpp"
#include "scalfmm/tree/group_tree.hpp"
#include "scalfmm/tree/leaf.hpp"
#include "scalfmm/utils/sort.hpp"
-
+#include "scalfmm/tools/colorized.hpp"
#include "scalfmm/tools/fma_loader.hpp"
#include "scalfmm/tools/vtk_writer.hpp"
#include "scalfmm/utils/accurater.hpp"
@@ -31,7 +26,9 @@
#include <chrono>
#include <cstdio>
#include <inria/tcli/tcli.hpp>
+#include <inria/tcli/help_descriptor.hpp>
#include <iostream>
+#include <limits>
#include <sstream>
#include <string>
#include <sys/types.h>
@@ -40,6 +37,13 @@
#include <unistd.h>
#include <utility>
#include <vector>
+
+//#define USE_CATALYST
+
+#ifdef USE_CATALYST
+#include "catalyst_related.hpp"
+#endif
+
// example:
// ./examples/RelWithDebInfo/test-time-loop --tree-height 4 -gs 2 --order 4 --dimension 2...
namespace local_args
@@ -90,7 +94,7 @@ namespace local_args
using type = int;
type def = 250;
};
- struct dimension
+ struct dimension : inria::tcli::required_tag
{
inria::tcli::str_vec flags = {"--dimension", "--d"};
const char* description = "Dimension : \n 2 for dimension 2, 3 for dimension 3";
@@ -99,13 +103,28 @@ namespace local_args
};
struct time_steps
{
- inria::tcli::str_vec flags = {"--time-step", "--ts"};
+ inria::tcli::str_vec flags = {"--time-steps", "--ts"};
const char* description = "Number of time steps";
using type = std::size_t;
type def = 1;
};
+ struct delta
+ {
+ inria::tcli::str_vec flags = {"--delta"};
+ const char* description = "Delta to apply on forces";
+ using type = double;
+ type def = 0.1;
+ };
+ struct catalyst
+ {
+ inria::tcli::str_vec flags = {"--catalyst"};
+ const char* description = "Enable catalist";
+ using type = bool;
+ type def = false;
+ };
} // 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)
@@ -141,6 +160,58 @@ void read_data(const std::string& filename, ContainerType*& container,
}
}
+template<typename GroupTree, typename ValueType>
+auto update_particles(GroupTree& tree, ValueType delta) -> void
+{
+ scalfmm::component::for_each_leaf(tree.begin(), tree.end(), [delta](auto& leaf)
+ {
+ auto& particles{leaf.particles()};
+ scalfmm::container::point<ValueType, GroupTree::dimension> force{};
+
+ for(std::size_t i{0}; i<particles.size(); ++i)
+ {
+ auto proxy = particles.at(i);
+ force = -delta*proxy.position();
+ proxy.position() += force;
+ }
+ });
+}
+
+template<typename GroupTree>
+auto get_new_box(GroupTree const& tree) -> typename GroupTree::box_type
+{
+ using value_type = typename GroupTree::leaf_type::particle_type::position_type::value_type;
+ using box_type = typename GroupTree::box_type;
+ auto const old_center = tree.box().center();
+
+ std::vector<value_type> max(GroupTree::dimension, 0);
+ std::vector<value_type> min(GroupTree::dimension, 0);
+
+ scalfmm::component::for_each_leaf(tree.begin(), tree.end(), [&min, &max, &old_center](auto const& leaf)
+ {
+ auto& particles{leaf.particles()};
+ for(std::size_t i{0}; i<particles.size(); ++i)
+ {
+ const auto position = particles.at(i).position();
+ //std::cout << particles.at(i).inputs(0) << '\n';
+ for(std::size_t d{0}; d<GroupTree::dimension; ++d)
+ {
+ min.at(d) = std::min(position.at(d), min.at(d));
+ max.at(d) = std::max(position.at(d), max.at(d));
+ //std::cout << particles.at(i).outputs(d) << ' ';
+ }
+ //std::cout << '\n';
+ }
+ });
+ value_type width{0};
+ for(std::size_t i{0}; i < GroupTree::dimension; ++i)
+ {
+ width = std::max(std::abs(min.at(i)), max.at(i));
+ }
+
+ return box_type(width + std::numeric_limits<value_type>::epsilon(), old_center);
+}
+
template<std::size_t Dimension, typename FmmOperatorType, typename... Parameters>
auto run(inria::tcli::parser<Parameters...> const& parser) -> int
{
@@ -148,6 +219,8 @@ auto run(inria::tcli::parser<Parameters...> const& parser) -> int
// timer
scalfmm::utils::timer time{};
+ static constexpr auto dimension{Dimension};
+
using value_type = double;
// ---------------------------------------
using near_matrix_kernel_type = typename FmmOperatorType::near_field_type::matrix_kernel_type;
@@ -179,11 +252,14 @@ auto run(inria::tcli::parser<Parameters...> const& parser) -> int
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>());
+ const bool catalyst_enable(parser.template get<local_args::catalyst>());
+ const auto delta = parser.template get<local_args::delta>();
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';
+ std::cout << scalfmm::colors::blue << "<params> Delta forces : " << delta << scalfmm::colors::reset << '\n';
// ---------------------------------------
// Open particle file
@@ -232,23 +308,45 @@ auto run(inria::tcli::parser<Parameters...> const& parser) -> int
std::cout << scalfmm::colors::yellow << "Kernel and Interp created in " << time.elapsed() << "ms\n"
<< scalfmm::colors::reset;
+#ifdef USE_CATALYST
+ if(catalyst_enable)
+ {
+ std::cout << scalfmm::colors::on_blue << "Initializing catalyst..." << scalfmm::colors::reset << '\n';
+ catalyst_adaptor::initialize(parser);
+ }
+#endif
// ---------------------------------------
// start algorithm for step 0
auto operator_to_proceed = scalfmm::algorithms::all;
- scalfmm::algorithms::sequential(tree, std::move(fmm_operator), operator_to_proceed);
+ scalfmm::algorithms::sequential(tree, fmm_operator, operator_to_proceed);
// ---------------------------------------
// time loop
for(std::size_t steps{1}; steps < n_time_steps; ++steps)
{
- // ---
- // compute the forces
// ---
// update particles
+ update_particles(tree, delta);
+ // ---
+ // visu related post treatment
+ if(!visu_file.empty())
+ {
+ scalfmm::tools::io::exportVTKxml(std::to_string(steps-1)+visu_file, tree, container->size());
+ }
+
+#ifdef USE_CATALYST
+ if(catalyst_enable)
+ {
+ catalyst_adaptor::execute(steps-1, tree);
+ }
+#endif
// ---
// compute the corners and set new box
+ auto const new_box = get_new_box(tree);
+ std::cout << scalfmm::colors::on_green << "New box width is " << new_box.width(0) << scalfmm::colors::reset << '\n';
// ---
// compute new morton indices
+
// ---
// move particles according to new morton indices
// ---
@@ -256,8 +354,29 @@ auto run(inria::tcli::parser<Parameters...> const& parser) -> int
// ---
// build tree levels
// ---
+ // start algorithm at step : compute forces
+ // ---
}
+ // ---------------------------------------
+ // last step update
+ // ---
+ // update particles
+ // ---
+ update_particles(tree, delta);
+
+ // ---
+ // visu related post treatment
+ if(!visu_file.empty())
+ {
+ scalfmm::tools::io::exportVTKxml(std::to_string(n_time_steps-1) + visu_file, tree, container->size());
+ }
+#ifdef USE_CATALYST
+ if(catalyst_enable)
+ {
+ }
+#endif
+
// ---------------------------------------
// write output fma
if(!output_file.empty())
@@ -283,10 +402,10 @@ 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{});
+ 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{}, local_args::delta{}, local_args::catalyst{});
parser.parse(argc, argv);
// Getting command line parameters
@@ -297,28 +416,29 @@ auto main([[maybe_unused]] int argc, [[maybe_unused]] char* argv[]) -> int
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 near_matrix_kernel_type = scalfmm::matrix_kernels::others::grad_one_over_r2<dim>;
+ using near_field_type = scalfmm::operators::near_field_operator<near_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>;
+ using far_matrix_kernel_type = scalfmm::matrix_kernels::others::one_over_r2;
+ using interpolator_type = scalfmm::interpolation::uniform_interpolator<double, dim, far_matrix_kernel_type>;
+ using far_field_type = scalfmm::operators::far_field_operator<interpolator_type, true>;
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>;
+ //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;
- }
+ // 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 "