diff --git a/bhtree_mpi/src/datastructures/BarnesHutTree.cpp b/bhtree_mpi/src/datastructures/BarnesHutTree.cpp index d1000a7e7f84190b4b7daf3709ebfa83f96d84fb..0ad4b48c09d04c71356e107f837238653c1efeae 100644 --- a/bhtree_mpi/src/datastructures/BarnesHutTree.cpp +++ b/bhtree_mpi/src/datastructures/BarnesHutTree.cpp @@ -23,33 +23,26 @@ namespace nbody { while (current != nodes) { current->representative.id = std::numeric_limits::max(); current->representative.mass = 0.0; - current->representative.position[0] = 0.0; - current->representative.position[1] = 0.0; - current->representative.position[2] = 0.0; + current->representative.position = Vec3{ 0.0 }; if (current->leaf) { for (std::size_t i = 0; i < current->bodies.size(); i++) { current->representative.mass += current->bodies[i].mass; - for (std::size_t j = 0; j < 3; j++) { - current->representative.position[j] += current->bodies[i].position[j] * current->bodies[i].mass; - } + current->representative.position += current->bodies[i].position * current->bodies[i].mass; } } else { Node* child = current->next; do { current->representative.mass += child->representative.mass; - for (std::size_t j = 0; j < 3; j++) { - current->representative.position[j] += child->representative.position[j] * child->representative.mass; - } + current->representative.position += child->representative.position * child->representative.mass; child = child->nextSibling; } while (child != nullptr); } - for (std::size_t j = 0; j < 3; j++) { - if (current->representative.mass > std::numeric_limits::epsilon()) { - current->representative.position[j] /= current->representative.mass; - } else { - current->representative.position[j] = 0.0; - } + if (current->representative.mass > std::numeric_limits::epsilon()) { + current->representative.position /= current->representative.mass; + } + else { + current->representative.position = Vec3{ 0.0 }; } current = current->prev; } @@ -61,19 +54,19 @@ namespace nbody { current->leaf = false; Node* after = current->next; - for (auto it = std::begin(subboxes); it != std::end(subboxes); it++) { + bool first = true; + for (const auto& subBox : subboxes) { Node* child = new Node(current->tree); - child->bb = *it; - child->bodies = it->copyBodies(current->bodies); - child->nextSibling = nullptr; - child->prevSibling = nullptr; + child->bb = subBox; + child->bodies = subBox.copyBodies(current->bodies); after->insertBefore(child); - if (it != std::begin(subboxes)) { + if (!first) { child->prev->nextSibling = child; child->prevSibling = child->prev; child->prev->afterSubtree = child; } + first = false; } after->prev->afterSubtree = current->afterSubtree; current->bodies.clear(); @@ -128,7 +121,7 @@ namespace nbody { //put all new bodies into fitting leaves, walk through tree and split Node* current; - for (auto& b : bodies) { + for (const auto& b : bodies) { current = nodes->next; while (!current->leaf) { Node* child = current->next; diff --git a/bhtree_mpi/src/datastructures/Body.cpp b/bhtree_mpi/src/datastructures/Body.cpp index 73a9e35baac9a61eeba9502761adfee43c2694cb..8ab72436121baa84a61a045d872c407ae46cf0fb 100644 --- a/bhtree_mpi/src/datastructures/Body.cpp +++ b/bhtree_mpi/src/datastructures/Body.cpp @@ -5,6 +5,7 @@ #include #include #include +#include namespace nbody { Body::Body(const Vec3& position_, @@ -34,41 +35,30 @@ namespace nbody { //integrating acceleration -> velocity -> position void Body::integrate() { - Derivative start; - Derivative a, b, c, d; - double dxdt[3]; - double dvdt[3]; - if (refinement) { return; } - a = evaluate(0.0, start); - b = evaluate(timestep * 0.5, a); - c = evaluate(timestep * 0.5, b); - d = evaluate(timestep, c); - for (std::size_t i = 0; i < 3; i++) { - dxdt[i] = 1.0 / 6.0 * (a.dx[i] + 2.0f * (b.dx[i] + c.dx[i]) + d.dx[i]); - dvdt[i] = 1.0 / 6.0 * (a.dv[i] + 2.0f * (b.dv[i] + c.dv[i]) + d.dv[i]); - position[i] = position[i] + dxdt[i] * timestep; - velocity[i] = velocity[i] + dvdt[i] * timestep; - } + const auto a = evaluate(0.0, Derivative{}); + const auto b = evaluate(timestep * 0.5, a); + const auto c = evaluate(timestep * 0.5, b); + const auto d = evaluate(timestep, c); + + const auto dxdt = Vec3{ 1.0 / 6.0 } *(a.dx + Vec3{ 2.0 } *(b.dx + c.dx) + d.dx); + const auto dvdt= Vec3{ 1.0 / 6.0 } * (a.dv + Vec3{ 2.0 } * (b.dv + c.dv) + d.dv); + position = position + dxdt * timestep; + velocity = velocity + dvdt * timestep; } //gravitational force computation void Body::accumulateForceOntoBody(const Body& from) { - double distance2 = 0.0; - - distance2 += (from.position[0] - position[0]) * (from.position[0] - position[0]); - distance2 += (from.position[1] - position[1]) * (from.position[1] - position[1]); - distance2 += (from.position[2] - position[2]) * (from.position[2] - position[2]); - distance2 = std::max(distance2, std::numeric_limits::epsilon()); - double distance = sqrt(distance2); - double mdist = -1.0 * ((from.mass * mass) / distance2); - for (std::size_t i = 0; i < 3; i++) { - double vec = (from.position[i] - position[i]) / distance; - - acceleration[i] += vec * mdist; - } + const auto distance2 = std::max( + SquaredDistance(from.position, position), + std::numeric_limits::epsilon()); + const auto distance = sqrt(distance2); + const auto mdist = -1.0 * ((from.mass * mass) / distance2); + + const auto vec = (from.position - position) / distance; + acceleration += vec * mdist; } bool validDouble(double val) { @@ -84,9 +74,9 @@ namespace nbody { } void Body::print(std::size_t parallelId) const { - std::cout << parallelId << " " << id << " Position: " << position[0] << " " << position[1] << " " << position[2] << '\n'; - std::cout << parallelId << " " << id << " Velocity: " << velocity[0] << " " << velocity[1] << " " << velocity[2] << '\n'; - std::cout << parallelId << " " << id << " Acceleration: " << acceleration[0] << " " << acceleration[1] << " " << acceleration[2] << '\n'; + std::cout << parallelId << " " << id << " Position: " << position << '\n'; + std::cout << parallelId << " " << id << " Velocity: " << velocity << '\n'; + std::cout << parallelId << " " << id << " Acceleration: " << acceleration << '\n'; std::cout << parallelId << " " << id << " Mass: " << mass << '\n'; std::cout << parallelId << " " << id << " Refinement: " << refinement << '\n' << '\n'; } diff --git a/bhtree_mpi/src/datastructures/Node.cpp b/bhtree_mpi/src/datastructures/Node.cpp index 35f1921ab1d434c61397515f7afdb6e5ae206aab..272763932200a8499db99d21ee8792f644243532 100644 --- a/bhtree_mpi/src/datastructures/Node.cpp +++ b/bhtree_mpi/src/datastructures/Node.cpp @@ -68,11 +68,9 @@ namespace nbody { std::cerr << "bb wrong\n"; return false; } - for (std::size_t i = 0; i < 3; i++) { - if (bb.min[i] > bb.max[i]) { - std::cerr << "bb " << i << " min " << bb.min[i] << " max " << bb.max[i] << '\n'; - return false; - } + if (bb.min > bb.max) { + std::cerr << "bb min " << bb.min << " max " << bb.max << '\n'; + return false; } if (std::any_of(std::begin(bodies), std::end(bodies), [&](const Body& b) {return !isContained(b, bb); })) { std::cerr << "bb out of bounds\n"; @@ -116,24 +114,20 @@ namespace nbody { //update overall node information void Node::update() { - double position[3] = {0.0, 0.0, 0.0}; + Vec3 position{0.0}; double mass = 0.0; if (leaf) { for (const auto& body : bodies) { mass += body.mass; - for (std::size_t i = 0; i < 3; i++) { - position[i] += body.position[i] * body.mass; - } + position += body.position * body.mass; } } else { for (Node* node = next; node != nullptr; node = node->nextSibling) { mass += node->representative.mass; } } - for (std::size_t i = 0; i < 3; i++) { - representative.position[i] = position[i] / mass; - } + representative.position = position / mass; representative.mass = mass; } @@ -152,12 +146,8 @@ namespace nbody { //check if node is sufficient for force evaluation bool Node::sufficientForBody(const Body& body) const { - double distance = 0.0; - - for (std::size_t i = 0; i < 3; i++) { - distance += (representative.position[i] - body.position[i]) * (representative.position[i] - body.position[i]); - } - return distance > std::pow(getL(), 2.0); + const auto distance2 = SquaredDistance(representative.position, body.position); + return distance2 > std::pow(getL(), 2.0); } //check if node is sufficient for force evaluation for all bodies in box diff --git a/bhtree_mpi/src/datastructures/Vec3.hpp b/bhtree_mpi/src/datastructures/Vec3.hpp index 8b582e4fdcce0fe1b93fecb3f28d3ff5025f6522..d7c20836d3ba3cf61fef75c6f21ffcaf161ca85f 100644 --- a/bhtree_mpi/src/datastructures/Vec3.hpp +++ b/bhtree_mpi/src/datastructures/Vec3.hpp @@ -156,9 +156,10 @@ inline bool operator==(const Vec3& lhs, const Vec3& rhs) { lhs.Z == rhs.Z; } -inline double Norm(const Vec3& v) { return std::sqrt(v.X * v.X + v.Y * v.Y + v.Z * v.Z); } - +inline double SquaredNorm(const Vec3& v) { return v.X * v.X + v.Y * v.Y + v.Z * v.Z; } +inline double Norm(const Vec3& v) { return std::sqrt(SquaredNorm(v)); } inline double Distance(const Vec3& v1, const Vec3& v2) { return Norm(v1 - v2); } +inline double SquaredDistance(const Vec3& v1, const Vec3& v2) { return SquaredNorm(v1 - v2); } inline std::ostream& operator<<(std::ostream& os, const Vec3 v) { os << v.X << ":" << v.Y << ":" << v.Z; diff --git a/bhtree_mpi/src/simulation/MpiSimulation.cpp b/bhtree_mpi/src/simulation/MpiSimulation.cpp index 16dfb0417ac2bf25d741745fa07048876c7c4199..27b1a611e29940ba35d5fab1a7ffdad68d100c1a 100644 --- a/bhtree_mpi/src/simulation/MpiSimulation.cpp +++ b/bhtree_mpi/src/simulation/MpiSimulation.cpp @@ -25,21 +25,23 @@ namespace nbody { int bodyBlocklengths[6] = { 1, 1, 1, 1, 1, 1 }; MPI_Datatype bodyDatatypes[6] = { MPI_UINT64_T, vec3Type, vec3Type, vec3Type, MPI_DOUBLE, MPI_CXX_BOOL }; - MPI_Aint bodyOffsets[6]; - bodyOffsets[0] = offsetof(Body, id); - bodyOffsets[1] = offsetof(Body, position); - bodyOffsets[2] = offsetof(Body, velocity); - bodyOffsets[3] = offsetof(Body, acceleration); - bodyOffsets[4] = offsetof(Body, mass); - bodyOffsets[5] = offsetof(Body, refinement); + MPI_Aint bodyOffsets[6] = { + offsetof(Body, id) , + offsetof(Body, position) , + offsetof(Body, velocity) , + offsetof(Body, acceleration), + offsetof(Body, mass) , + offsetof(Body, refinement) + }; MPI_Type_create_struct(6, bodyBlocklengths, bodyOffsets, bodyDatatypes, &bodyType); MPI_Type_commit(&bodyType); int boxBlocklengths[2] = { 1, 1 }; MPI_Datatype boxDatatypes[2] = { vec3Type, vec3Type }; - MPI_Aint boxOffsets[2]; - boxOffsets[0] = offsetof(Box, min); - boxOffsets[1] = offsetof(Box, max); + MPI_Aint boxOffsets[2] = { + offsetof(Box, min), + offsetof(Box, max) + }; MPI_Type_create_struct(2, boxBlocklengths, boxOffsets, boxDatatypes, &boxType); MPI_Type_commit(&boxType);