Skip to content

Commit

Permalink
Organize different gravity derivative implementations, including Chri…
Browse files Browse the repository at this point in the history 
…stoph's generated code.
  • Loading branch information
@JacksonCampolattaro
JacksonCampolattaro committed Jun 29, 2024
1 parent c85ba55 commit f1463b4
Show file tree
Hide file tree
Showing 9 changed files with 847 additions and 160 deletions.
13 changes: 6 additions & 7 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,15 +7,14 @@
![MIT license](https://img.shields.io/badge/license-MIT-A51931)
![C++20](https://img.shields.io/badge/standard-C%2B%2B20-blue?logo=cplusplus)

_For more information, see the [documentation](https://jackson.campolattaro.nl/symtensor)!_
_For more information, see the [documentation](https://jackson.campolattaro.nl/symtensor)._

Symtensor is designed for working with tensors which are...
- **Symmetric** — only the upper triangle needs to be stored.
- **Small** — indices have semantic meaning, similar to GLM's `mat3x3` type.
- **High-rank** — not just vectors and matrices, these tensors can have any number of dimensions.
Symtensor is designed for working with tensors which are **Symmetric**, **Small** and **High-rank**, and for the composition of those tensors into Multipoles.

These are unusual needs for a tensor library. The design comes from symtensor's original purpose as part of an implementation of the Fast Multipole Method for an [_n_-body gravity simulation](https://github.com/JacksonCampolattaro/n-body). Popular tensor libraries like [GLM](https://github.com/g-truc/glm), [Eigen](https://github.com/libigl/eigen), and [Fastor](https://github.com/romeric/Fastor) are immediately ruled out for this purpose.
These are unusual needs for a tensor library. The design comes from symtensor's original purpose as part of an implementation of the Fast Multipole Method for an [_n_-body gravity simulation](https://github.com/JacksonCampolattaro/n-body). Popular tensor libraries like [GLM](https://github.com/g-truc/glm), [Eigen](https://github.com/libigl/eigen), and [Fastor](https://github.com/romeric/Fastor) are ruled out for this purpose because they don't support condensed storage of symmetric tensors.

There are a handful of projects which aim to solve the same problems as this one, including [Tensor](https://github.com/thenumbernine/Tensor), [GADGET-4's symmetric tensor library](https://github.com/weiguangcui/Gadget4/blob/Gadget4-Simba/src/data/symtensors.h), and [fmmgen](https://github.com/rpep/fmmgen). Tensor is feature-rich, but has a broader maths focus and provides no benchmarks. GADGET-4's library has been proven to perform well for _n_-body simulations, but is ultimately limited by its handwritten operators -- it can't be used for high-order multipoles. Fmmgen uses a script to generate efficient multipole and symmetric tensor code which makes it very flexible, but the code it produces is entirely based on C-style arrays and not intended to be human-readable.
There are a handful of projects which aim to solve the same problems as this one, including [Tensor](https://github.com/thenumbernine/Tensor), [GADGET-4's symmetric tensor library](https://github.com/weiguangcui/Gadget4/blob/Gadget4-Simba/src/data/symtensors.h), and [fmmgen](https://github.com/rpep/fmmgen). Tensor is feature-rich, but has a broader maths focus and doesn't provide benchmarks. GADGET-4's library has been proven to perform well for _n_-body simulations, but is ultimately limited by its handwritten operators -- it can't be used for high-order multipoles. Fmmgen uses a script to generate efficient multipole and symmetric tensor code which makes it very flexible, but the code it produces is based on C-style arrays and not intended to be human-readable.

The goal of this project is to produce a library which combines the modern sensibilities of _Tensor_ with the performance of _GADGET-4_'s library and the flexibility of _fmmgen_.

Eventually, I intend to incorporate symtensor into a larger _n_-body solver built from the ground up. To this end, [Christoph Peters](https://momentsingraphics.de/About.html) has been a big help by generating derivatives of the Newton's gravity equation up to arbitrarily high orders.
5 changes: 5 additions & 0 deletions benchmarks/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,3 +27,8 @@ target_link_libraries(benchmarks PRIVATE
Catch2::Catch2WithMain
symtensor::symtensor
)
target_compile_options(
benchmarks PRIVATE
$<$<CXX_COMPILER_ID:Clang>:-fassociative-math>
$<$<CXX_COMPILER_ID:AppleClang>:-fassociative-math>
)
67 changes: 39 additions & 28 deletions benchmarks/multipoleMoment.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,10 @@
#include <glm/geometric.hpp>

#include <symtensor/Multipole.h>
#include <symtensor/gravity.h>

#include "symtensor/gravity/direct.h"
#include "symtensor/gravity/einsum.h"
#include "symtensor/gravity/tensorlib.h"

using namespace symtensor;

Expand All @@ -16,31 +19,39 @@ TEST_CASE("benchmark: Multipole moment gravity approximation", "[MultipoleMoment
TEST_CASE("benchmark: Gravity derivatives construction", "[Gravity]") {
auto R = glm::vec3{1.0, 2.0, 3.0};


constexpr std::array<int, 5> input = {1, 2, 2, 3, 1};
constexpr auto result = filter<input, [&](auto v) constexpr { return v != 3; }>();

// Print the result (for non-constexpr context)
for (const auto &value: result) {
std::cout << value << std::endl;
}

CHECK((gravity::D<2>(R) - gravity::derivative<2>(R)).norm() < 1e-7);
CHECK((gravity::D<3>(R) - gravity::derivative<3>(R)).norm() < 1e-7);
CHECK((gravity::D<4>(R) - gravity::derivative<4>(R)).norm() < 1e-7);
CHECK((gravity::D<5>(R) - gravity::derivative<5>(R)).norm() < 1e-7); // todo

// BENCHMARK("D' Construction") { return gravity::D<1>(R); };
// BENCHMARK("D' Construction (improved)") { return gravity::derivative<1>(R); };
// BENCHMARK("D'' Construction") { return gravity::D<2>(R); };
// BENCHMARK("D'' Construction (improved)") { return gravity::derivative<2>(R); };
// BENCHMARK("D''' Construction") { return gravity::D<3>(R); };
// BENCHMARK("D''' Construction (improved)") { return gravity::derivative<3>(R); };
// BENCHMARK("D'''' Construction") { return gravity::D<4>(R); };
// BENCHMARK("D'''' Construction (improved)") { return gravity::derivative<4>(R); };
// BENCHMARK("D''''' Construction") { return gravity::D<5>(R); };
// BENCHMARK("D''''' Construction (improved)") { return gravity::derivative<5>(R); }; // todo

BENCHMARK("D1-4 Construction") { return gravity::Ds<4>(R); };
BENCHMARK("D1-4 Construction (improved)") { return gravity::derivatives<4>(R); };
// Einsum must match direct form
CHECK((gravity::direct::derivative<2>(R) - gravity::einsum::derivative<2>(R)).norm() < 1e-7);
CHECK((gravity::direct::derivative<3>(R) - gravity::einsum::derivative<3>(R)).norm() < 1e-7);
CHECK((gravity::direct::derivative<4>(R) - gravity::einsum::derivative<4>(R)).norm() < 1e-7);
// CHECK((gravity::D<5>(R) - gravity::derivative<5>(R)).norm() < 1e-7); // todo: unimplemented

// Tensorlib must match direct form
CHECK((gravity::direct::derivative<2>(R) - gravity::tensorlib::derivative<2>(R)).norm() < 1e-7);
CHECK((gravity::direct::derivative<3>(R) - gravity::tensorlib::derivative<3>(R)).norm() < 1e-7);
CHECK((gravity::direct::derivative<4>(R) - gravity::tensorlib::derivative<4>(R)).norm() < 1e-7);
CHECK((gravity::direct::derivative<5>(R) - gravity::tensorlib::derivative<5>(R)).norm() < 1e-7);

BENCHMARK("D' (direct)") { return gravity::direct::derivative<1>(R); };
BENCHMARK("D' (einsum)") { return gravity::einsum::derivative<1>(R); };
BENCHMARK("D' (tensorlib)") { return gravity::tensorlib::derivative<1>(R); };
BENCHMARK("D'' (direct)") { return gravity::direct::derivative<2>(R); };
BENCHMARK("D'' (einsum)") { return gravity::einsum::derivative<2>(R); };
BENCHMARK("D'' (tensorlib)") { return gravity::tensorlib::derivative<2>(R); };
BENCHMARK("D''' (direct)") { return gravity::direct::derivative<3>(R); };
BENCHMARK("D''' (einsum)") { return gravity::einsum::derivative<3>(R); };
BENCHMARK("D''' (tensorlib)") { return gravity::tensorlib::derivative<3>(R); };
BENCHMARK("D'''' (direct)") { return gravity::direct::derivative<4>(R); };
BENCHMARK("D'''' (einsum)") { return gravity::einsum::derivative<4>(R); };
BENCHMARK("D'''' (tensorlib)") { return gravity::tensorlib::derivative<4>(R); };

BENCHMARK("D' - D'''' (direct)") { return gravity::direct::derivatives<4>(R); };
BENCHMARK("D' - D'''' (einsum)") { return gravity::einsum::derivatives<4>(R); };
BENCHMARK("D' - D'''' (tensorlib)") { return gravity::tensorlib::derivatives<4>(R); };

BENCHMARK("D' - D''''' (direct)") { return gravity::direct::derivatives<5>(R); };
//BENCHMARK("D' - D''''' (einsum)") { return gravity::einsum::derivatives<5>(R); }; // todo: unimplemented
BENCHMARK("D' - D''''' (tensorlib)") { return gravity::tensorlib::derivatives<5>(R); };

// BENCHMARK("D1-5 Construction (direct)") { return gravity::direct::derivatives<5>(R); };
// BENCHMARK("D1-5 Construction (tensorlib)") { return gravity::tensorlib::derivatives<5>(R); };
}
13 changes: 6 additions & 7 deletions doc/pages/main.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,13 +11,12 @@
![C++20](https://img.shields.io/badge/standard-C%2B%2B20-blue?logo=cplusplus)
@m_enddiv

Symtensor is designed for working with tensors which are...
- **Symmetric** — only the upper triangle needs to be stored.
- **Small** — indices have semantic meaning, similar to GLM's `mat3x3` type.
- **High-rank** — not just vectors and matrices, these tensors can have any number of dimensions.
Symtensor is designed for working with tensors which are **Symmetric**, **Small** and **High-rank**, and for the composition of those tensors into Multipoles.

These are unusual needs for a tensor library. The design comes from symtensor's original purpose as part of an implementation of the Fast Multipole Method for an [_n_-body gravity simulation](https://github.com/JacksonCampolattaro/n-body). Popular tensor libraries like [GLM](https://github.com/g-truc/glm), [Eigen](https://github.com/libigl/eigen), and [Fastor](https://github.com/romeric/Fastor) are immediately ruled out for this purpose.
These are unusual needs for a tensor library. The design comes from symtensor's original purpose as part of an implementation of the Fast Multipole Method for an [_n_-body gravity simulation](https://github.com/JacksonCampolattaro/n-body). Popular tensor libraries like [GLM](https://github.com/g-truc/glm), [Eigen](https://github.com/libigl/eigen), and [Fastor](https://github.com/romeric/Fastor) are ruled out for this purpose because they don't support condensed storage of symmetric tensors.

There are a handful of projects which aim to solve the same problems as this one, including [Tensor](https://github.com/thenumbernine/Tensor), [GADGET-4's symmetric tensor library](https://github.com/weiguangcui/Gadget4/blob/Gadget4-Simba/src/data/symtensors.h), and [fmmgen](https://github.com/rpep/fmmgen). Tensor is feature-rich, but has a broader maths focus and provides no benchmarks. GADGET-4's library has been proven to perform well for _n_-body simulations, but is ultimately limited by its handwritten operators -- it can't be used for high-order multipoles. Fmmgen uses a script to generate efficient multipole and symmetric tensor code which makes it very flexible, but the code it produces is entirely based on C-style arrays and not intended to be human-readable.
There are a handful of projects which aim to solve the same problems as this one, including [Tensor](https://github.com/thenumbernine/Tensor), [GADGET-4's symmetric tensor library](https://github.com/weiguangcui/Gadget4/blob/Gadget4-Simba/src/data/symtensors.h), and [fmmgen](https://github.com/rpep/fmmgen). Tensor is feature-rich, but has a broader maths focus and doesn't provide benchmarks. GADGET-4's library has been proven to perform well for _n_-body simulations, but is ultimately limited by its handwritten operators -- it can't be used for high-order multipoles. Fmmgen uses a script to generate efficient multipole and symmetric tensor code which makes it very flexible, but the code it produces is based on C-style arrays and not intended to be human-readable.

The goal of this project is to produce a library which combines the modern sensibilities of _Tensor_ with the performance of _GADGET-4_'s library and the flexibility of _fmmgen_.
The goal of this project is to produce a library which combines the modern sensibilities of _Tensor_ with the performance of _GADGET-4_'s library and the flexibility of _fmmgen_.

Eventually, I intend to incorporate symtensor into a larger _n_-body solver built from the ground up. To this end, [Christoph Peters](https://momentsingraphics.de/About.html) has been a big help by generating derivatives of the Newton's gravity equation up to arbitrarily high orders.
25 changes: 16 additions & 9 deletions include/symtensor/SymmetricTensorBase.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -167,6 +167,13 @@ namespace symtensor {
static_assert(sizeof...(s) == NumUniqueValues);
}

/**
* @brief Constructor from an std::array of scalars.
*
* @param s an array of scalar values to initialize the tensor.
*/
explicit constexpr SymmetricTensorBase(const std::array<S, NumUniqueValues> &values) : _data{values} {}

/**
* @brief (Implicit) constructor from a sequence of values.
*
Expand Down Expand Up @@ -692,15 +699,15 @@ namespace symtensor {
static_assert(NumValues > 0);

// Use a lookup table when the range of options is small enough
// if constexpr (R < 3)
// return [&]<std::size_t... i>(std::index_sequence<i...>) {
// return as_lookup_table<
// decltype([](std::array<I, R> ind) consteval { return symtensor::flatIndex(ind, D); }),
// std::array<I, R>,
// lexicographicalIndices(i)...
// >(indices);
// }(std::make_index_sequence<NumValues>());
// else
// if constexpr (R < 3)
// return [&]<std::size_t... i>(std::index_sequence<i...>) {
// return as_lookup_table<
// decltype([](std::array<I, R> ind) consteval { return symtensor::flatIndex(ind, D); }),
// std::array<I, R>,
// lexicographicalIndices(i)...
// >(indices);
// }(std::make_index_sequence<NumValues>());
// else
return symtensor::flatIndex(indices, D);
}

Expand Down
120 changes: 11 additions & 109 deletions include/symtensor/gravity.h → include/symtensor/gravity/direct.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,104 +4,14 @@

#include <glm/gtx/norm.hpp>

#include <symtensor/symtensor.h>
#include <symtensor/util.h>
#include "symtensor/symtensor.h"
#include "symtensor/util.h"

namespace symtensor::gravity {
namespace symtensor::gravity::direct {

template<auto index, indexable Vector>
inline constexpr auto product_of_elements(const Vector &v) {
// todo: maybe don't use std::apply for this?
return std::apply([&](auto ...i) {
return (v[static_cast<std::size_t>(i)] * ...);
}, index);
}

template<auto index, std::size_t N, indexable Vector>
ALWAYS_INLINE auto derivative_at(const Vector &R) {

auto r2 = glm::length2(R);
auto r = sqrt(r2);
auto inv_r2 = 1.0f / r2;
auto g1 = -inv_r2 / r;
auto g2 = -3.0f * g1 * inv_r2;
auto g3 = -5.0f * g2 * inv_r2;
auto g4 = -7.0f * g3 * inv_r2;
auto g5 = -9.0f * g4 * inv_r2;

// auto r = glm::length(R);
// auto g1 = -1.0f / pow<3>(r);
// auto g2 = 3.0f / pow<5>(r);
// auto g3 = -15.0f / pow<7>(r);
// auto g4 = 105.0f / pow<9>(r);
// auto g5 = -945.0f / pow<11>(r);

auto cartesian_product_term = product_of_elements<index>(R);
if constexpr (N == 1) {
return R[static_cast<std::size_t>(index[0])] * g1;
} else if constexpr (N == 2) {
return g1 * kronecker_delta<float>(index) + g2 * cartesian_product_term;
} else if constexpr (N == 3) {
constexpr auto partitions = binomial_partitions<1>(index);
constexpr auto r_indices = std::get<0>(unzip(partitions));
constexpr auto kronecker_indices = std::get<1>(unzip(partitions));
constexpr auto repeats = repeats_table(kronecker_indices);
auto kronecker_product_term = [&]<auto... i>(std::index_sequence<i...>) LAMBDA_ALWAYS_INLINE {
return ((
force_consteval<kronecker_delta<kronecker_indices[i]>() && (repeats[i] > 0)> ?
product_of_elements<r_indices[i]>(R) * repeats[i] : 0
) + ...);
}(std::make_index_sequence<kronecker_indices.size()>());
return g2 * kronecker_product_term + g3 * cartesian_product_term;
} else if constexpr (N == 4) {
auto result = g4 * cartesian_product_term;
constexpr auto partitions = binomial_partitions<2>(index);
constexpr auto kronecker_indices = std::get<0>(unzip(partitions));
constexpr auto r_indices = std::get<1>(unzip(partitions));
{
constexpr auto kronecker_term = [&]<auto... i>(std::index_sequence<i...>) constexpr {
return ((kronecker_delta(kronecker_indices[i]) & kronecker_delta(r_indices[i])) + ...);
}(std::make_index_sequence<kronecker_indices.size() / 2>());
result += g2 * kronecker_term;
}
{
constexpr auto repeats = repeats_table(kronecker_indices);
auto kronecker_product_term = [&]<auto... i>(std::index_sequence<i...>) LAMBDA_ALWAYS_INLINE {
return ((
force_consteval<kronecker_delta<kronecker_indices[i]>() && (repeats[i] > 0)> ?
product_of_elements<r_indices[i]>(R) * repeats[i] : 0
) + ...);
}(std::make_index_sequence<kronecker_indices.size()>());
result += g3 * kronecker_product_term;
}

return result;
} else if constexpr (N == 5) {
auto result = g5 * cartesian_product_term;
return result;
}
return 0.0f;
}

template<std::size_t N, indexable Vector>
ALWAYS_INLINE auto derivative(const Vector &R) {
return SymmetricTensor3f<N>::NullaryExpression([&]<auto index>() LAMBDA_ALWAYS_INLINE {
return derivative_at<index, N>(R);
});
};

template<std::size_t N, indexable Vector>
ALWAYS_INLINE auto derivatives(const Vector &R) {
return [&]<auto... i>(std::index_sequence<i...>) LAMBDA_ALWAYS_INLINE {
// return Multipole3f<N>(derivative<i + 1>(R) ...);
return std::make_tuple(derivative<i + 1>(R) ...);
}(std::make_index_sequence<N>());
};


// todo: this needs to be regularized -- there must be a pattern!
// An implementation of gravity derivatives in the style of GADGET-4
template<std::size_t Order, indexable Vector>
ALWAYS_INLINE auto D(const Vector &R) {
ALWAYS_INLINE auto derivative(const Vector &R) {

auto r = glm::length(R);

Expand Down Expand Up @@ -181,10 +91,10 @@ namespace symtensor::gravity {
} else if constexpr (Order == 5) {


SymmetricTensor3f<2> R2 = SymmetricTensor3f<2>::CartesianPower(R);
SymmetricTensor3f<3> R3 = SymmetricTensor3f<3>::CartesianPower(R);
SymmetricTensor3f<2> R2 = SymmetricTensor3f < 2 > ::CartesianPower(R);
SymmetricTensor3f<3> R3 = SymmetricTensor3f < 3 > ::CartesianPower(R);

auto A = g5 * SymmetricTensor3f<5>::CartesianPower(R);
auto A = g5 * SymmetricTensor3f < 5 > ::CartesianPower(R);

// ~~~

Expand Down Expand Up @@ -251,19 +161,11 @@ namespace symtensor::gravity {
}

template<std::size_t N, indexable Vector>
ALWAYS_INLINE auto Ds(const Vector &R) {
ALWAYS_INLINE auto derivatives(const Vector &R) {
return [&]<auto... i>(std::index_sequence<i...>) LAMBDA_ALWAYS_INLINE {
return std::make_tuple(D<i + 1>(R) ...);
return std::make_tuple(derivative<i + 1>(R) ...);
}(std::make_index_sequence<N>());
};


auto derivative4(const glm::vec3 &R) {
return derivatives<4>(R);
}

auto D4(const glm::vec3 &R) {
return Ds<4>(R);
}

}
}
Loading

0 comments on commit f1463b4

Please sign in to comment.