Skip to content

Instantly share code, notes, and snippets.

@bernhardmgruber
Last active April 12, 2021 10:28
Show Gist options
  • Save bernhardmgruber/6a7a8d6b5fc65a0e1a566ba763f1bb6d to your computer and use it in GitHub Desktop.
Save bernhardmgruber/6a7a8d6b5fc65a0e1a566ba763f1bb6d to your computer and use it in GitHub Desktop.
n-body simulation using an SoAContainer using the proposed C++ reflection/injection facilities
#include <vector>
#include <tuple>
#include <cmath>
#include <random>
namespace meta = std::experimental::meta;
namespace lib {
namespace internal {
consteval void genSoRMembers(meta::info sov) {
for (meta::info member : meta::data_member_range(sov)) {
auto type = meta::type_of(member);
auto name = meta::name_of(member);
-> fragment struct {
typename(%{type})& unqualid(%{name});
};
}
}
template<typename T, typename U>
void assignEquallyNamedMembers(T& t, const U& u) {
consteval {
for (meta::info member : meta::data_member_range(reflexpr(T))) {
auto name = meta::name_of(member);
-> fragment {
idexpr(%{reflexpr(t)}).unqualid(%{name}) = idexpr(%{reflexpr(u)}).unqualid(%{name});
};
}
}
}
template <typename T>
struct SoR {
consteval {
-> fragment struct {
using type = struct Type {
consteval {
genSoRMembers(reflexpr(T));
}
operator T() const {
T t;
assignEquallyNamedMembers(t, *this);
return t;
}
auto operator=(const T& t) -> Type& {
assignEquallyNamedMembers(*this, t);
return *this;
}
};
};
}
};
template <typename T>
inline constexpr std::size_t memberCount = []() consteval {
auto i = 0;
for (meta::info member : meta::data_member_range(reflexpr(T)))
i++;
return i;
}();
template <typename T, std::size_t I>
struct MemberTypeAt {
consteval {
auto i = 0;
for (meta::info member : meta::data_member_range(reflexpr(T))) {
if (i == I) {
auto type = meta::type_of(member);
-> fragment struct {
using type = typename(%{type});
};
}
i++;
}
}
};
template <typename T, std::size_t... Is>
auto tupleOfVectorsType(std::index_sequence<Is...>) {
return std::tuple{(std::vector<typename MemberTypeAt<T, Is>::type>{})...};
}
}
template <typename T>
struct IndirectValue {
T value;
auto operator->() -> T* {
return &value;
}
auto operator->() const -> const T* {
return &value;
}
};
template <typename T>
struct SoAContainer {
using value_type = T;
using reference = typename internal::SoR<T>::type;
SoAContainer(std::size_t size) {
resize(size);
}
void resize(std::size_t size) {
resizeVectors(size, TupleIndexSeq{});
}
auto size() const -> std::size_t {
return std::get<0>(tupleOfVectors).size();
}
auto operator[](std::size_t i) -> reference {
return constructReferenceAt(i, TupleIndexSeq{});
}
auto operator[](std::size_t i) const -> value_type {
return constructValueAt(i, TupleIndexSeq{});
}
struct Iterator {
using iterator_category = std::random_access_iterator_tag;
using value_type = SoAContainer::value_type;
using difference_type = std::ptrdiff_t;
using pointer = IndirectValue<value_type>;
using reference = SoAContainer::reference;
constexpr auto operator++() -> Iterator& {
++pos;
return *this;
}
constexpr auto operator++(int) -> Iterator {
auto tmp = *this;
++*this;
return tmp;
}
constexpr auto operator--() -> Iterator& {
--pos;
return *this;
}
constexpr auto operator--(int) -> Iterator {
auto tmp{*this};
--*this;
return tmp;
}
constexpr auto operator*() const -> reference {
return (*cont)[pos];
}
constexpr auto operator->() const -> pointer {
return IndirectValue{**this};
}
constexpr auto operator[](difference_type i) const -> reference {
return *(*this + i);
}
constexpr auto operator+=(difference_type n) -> Iterator& {
pos = static_cast<difference_type>(pos) + n;
return *this;
}
friend constexpr auto operator+(Iterator it, difference_type n) -> Iterator {
it += n;
return it;
}
friend constexpr auto operator+(difference_type n, Iterator it) -> Iterator {
return it + n;
}
constexpr auto operator-=(difference_type n) -> Iterator& {
pos = static_cast<difference_type>(pos) - n;
return *this;
}
friend constexpr auto operator-(Iterator it, difference_type n) -> Iterator {
it -= n;
return it;
}
friend constexpr auto operator-(const Iterator& a, const Iterator& b) -> difference_type {
return static_cast<std::ptrdiff_t>(a.pos) - static_cast<std::ptrdiff_t>(b.pos);
}
friend constexpr auto operator==(const Iterator& a, const Iterator& b) -> bool {
return a.pos == b.pos;
}
friend constexpr auto operator!=(const Iterator& a, const Iterator& b) -> bool {
return !(a == b);
}
friend constexpr auto operator<(const Iterator& a, const Iterator& b) -> bool {
return a.pos < b.pos;
}
friend constexpr auto operator>(const Iterator& a, const Iterator& b) -> bool {
return b < a;
}
friend constexpr auto operator<=(const Iterator& a, const Iterator& b) -> bool {
return !(a > b);
}
friend constexpr auto operator>=(const Iterator& a, const Iterator& b) -> bool {
return !(a < b);
}
SoAContainer* cont;
std::size_t pos;
};
auto begin() -> Iterator {
return {this, 0};
}
auto end() -> Iterator {
return {this, size()};
}
private:
using TupleIndexSeq = std::make_index_sequence<internal::memberCount<T>>;
using TupleOfVectors = decltype(internal::tupleOfVectorsType<T>(TupleIndexSeq{}));
template <std::size_t... Is>
auto constructReferenceAt(std::size_t i, std::index_sequence<Is...>) {
return reference{std::get<Is>(tupleOfVectors)[i]...};
}
template <std::size_t... Is>
auto constructValueAt(std::size_t i, std::index_sequence<Is...>) const {
return value_type{std::get<Is>(tupleOfVectors)[i]...};
}
template <std::size_t... Is>
void resizeVectors(std::size_t size, std::index_sequence<Is...>) {
((std::get<Is>(tupleOfVectors).resize(size)),...);
}
TupleOfVectors tupleOfVectors;
};
}
constexpr auto PROBLEM_SIZE = 256;
constexpr auto TIMESTEP = 0.0001f;
constexpr auto EPS2 = 0.01f;
struct Particle {
float posx, posy, posz;
float velx, vely, velz;
float mass;
};
inline void pPInteraction(auto& pi, const auto& pj) {
auto xdistance = pi.posx - pj.posx;
auto ydistance = pi.posy - pj.posy;
auto zdistance = pi.posz - pj.posz;
xdistance *= xdistance;
ydistance *= ydistance;
zdistance *= zdistance;
const float distSqr = EPS2 + xdistance + ydistance + zdistance;
const float distSixth = distSqr * distSqr * distSqr;
const float invDistCube = 1.0f / std::sqrt(distSixth);
const float sts = pj.mass * invDistCube * TIMESTEP;
pi.velx += xdistance * sts;
pi.vely += ydistance * sts;
pi.velz += zdistance * sts;
}
void update(lib::SoAContainer<Particle>& particles) {
#pragma clang loop vectorize(enable)
for (std::size_t i = 0; i < PROBLEM_SIZE; i++) {
Particle pi = particles[i]; // SoR<Particle>::operator Particle()
//#pragma clang loop vectorize(enable)
for (auto&& pj : particles)
pPInteraction(pi, pj);
particles[i] = pi; // SoR<Particle>::operator=(Particle)
}
}
int main() {
lib::SoAContainer<Particle> c(PROBLEM_SIZE);
std::default_random_engine engine;
std::normal_distribution<float> dist(0.0f, 1.0f);
for (auto p : c) { // p is a reference
p.posx = dist(engine);
p.posy = dist(engine);
p.posz = dist(engine);
p.velx = dist(engine);
p.vely = dist(engine);
p.velz = dist(engine);
p.mass = dist(engine);
}
update(c);
return c[42].posz;
}
@bernhardmgruber
Copy link
Author

bernhardmgruber commented Mar 25, 2021

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment