Skip to content
Snippets Groups Projects
Commit f843e07d authored by Armin Damon Riess's avatar Armin Damon Riess
Browse files

remove member Q_, implement multipole (incorrect)

parent 34b555e8
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
lib/node.cpp
+ 96
− 71
View file @ f843e07d
......@@ -56,27 +56,6 @@ void Node::update(List localParticles, unsigned nLocalParticles) {
updateMassAndCOM();
// update matrix Q for multipole expansion
// go through every entry in Q
for (unsigned i=0; i<3; ++i) {
for (unsigned j=0; j<3; ++j) {
Q_[3*i+j] = 0;
// go through every local particle
for (unsigned c=0; c<8; ++c) {
if (children_[c] == nullptr) continue;
for (auto k : childParticles_[c]) {
Q_[3*i+j] += mass_ * 3*(centerOfMass_[i] - positions_[3*k+i])*(centerOfMass_[j] - positions_[3*k+j]);
if (i!=j) continue;
for (unsigned l=0; l<3; ++l) {
Q_[3*i+j] -= mass_ * (centerOfMass_[l] - positions_[3*k+l])*(centerOfMass_[l] - positions_[3*k+l]);
}
}
}
}
}
isLeaf_ = true;
if (depth_ >= maxDepth_) {
// if the node contains less than minNParticles, it is a leaf
......@@ -90,7 +69,7 @@ void Node::update(List localParticles, unsigned nLocalParticles) {
for (unsigned i=0; i<8; ++i) {
if (children_[i] != nullptr) {
// if child is not a nullptr and particleDivision[i] is empty or contains a single particle, delete the child
// if child is not a nullptr and particleDivision[i] contains too few particles, delete the child
// otherwise update the child
if (childParticles_[i].size() <= minNParticles_) {
delete children_[i];
......@@ -100,17 +79,18 @@ void Node::update(List localParticles, unsigned nLocalParticles) {
isLeaf_ = false;
}
} else {
// if child is a nullptr and particleDivision[i] contains more than one particle, create a new child
if (childParticles_[i].size() > minNParticles_) {
double childSize = size_/2;
Vector childCenter(3);
for (unsigned j=0; j<3; ++j) {
childCenter[j] = center_[j] + (i & (1 << j) ? childSize : -childSize);
}
children_[i] = new Node(mass_, positions_, nParticles_, softening_, depth_+1,
childSize, childCenter, childParticles_[i], childParticles_[i].size());
isLeaf_ = false;
// if child wouldn't contain enough particles, don't create a new child
if (childParticles_[i].size() < minNParticles_) continue;
// if child is a nullptr and particleDivision[i] contains enough particles, create a new child
double childSize = size_/2;
Vector childCenter(3);
for (unsigned j=0; j<3; ++j) {
childCenter[j] = center_[j] + (i & (1 << j) ? childSize : -childSize);
}
children_[i] = new Node(mass_, positions_, nParticles_, softening_, depth_+1,
childSize, childCenter, childParticles_[i], childParticles_[i].size());
isLeaf_ = false;
}
}
}
......@@ -124,60 +104,105 @@ Vector Node::calculateForce(unsigned particle) {
double dx = centerOfMass_[0] - px;
double dy = centerOfMass_[1] - py;
double dz = centerOfMass_[2] - pz;
double lambda = std::sqrt(dx*dx + dy*dy + dz*dz);
double theta = size_/lambda;
Vector y = {dx, dy, dz};
double ynorm = std::sqrt(dx*dx + dy*dy + dz*dz);
double theta = size_/ynorm;
// If theta small enough, calculate force on particle from this node with multipole expansion
if (theta < theta0_) {
force = multipole(y, ynorm);
return force;
}
// If theta too large:
// For each child that exists, calculate force on particle from that child.
// If child doesn't exist, calculate force on particle from particles in that child.
// If theta small enough, calculate force on particle from this node.
if (theta > theta0_) {
// get force for each child that exists and add them together
for (unsigned c=0; c<8; ++c) {
// if the child exists, get the force from the child
if (children_[c] != nullptr) {
Vector childForce = children_[c]->calculateForce(particle);
for (unsigned j=0; j<3; ++j) {
force[j] += childForce[j];
}
continue;
// If child doesn't exist, calculate force on particle from particles in that octant.
// get force for each child that exists and add them together
for (unsigned c=0; c<8; ++c) {
// if the child exists, get the force from the child
if (children_[c] != nullptr) {
Vector childForce = children_[c]->calculateForce(particle);
for (unsigned j=0; j<3; ++j) {
force[j] += childForce[j];
}
continue;
}
// if the child is a nullptr, calculate the force from the particles in the child (octant) directly
for (unsigned p : childParticles_[c]) {
// don't calculate force on particle from itself
if (p == particle) continue;
double pdx = positions_[3*p + 0] - px;
double pdy = positions_[3*p + 1] - py;
double pdz = positions_[3*p + 2] - pz;
double r2 = pdx*pdx + pdy*pdy + pdz*pdz;
double r = std::sqrt(r2);
// if the child is a nullptr, calculate the force from the particles in the child (octant) directly
for (unsigned p : childParticles_[c]) {
// don't calculate force on particle from itself
if (p == particle) continue;
double factor = mass_*mass_ / (r2 + softening_*softening_);
double pdx = positions_[3*p + 0] - px;
double pdy = positions_[3*p + 1] - py;
double pdz = positions_[3*p + 2] - pz;
// add the force from each particle in the child
force[0] += factor * pdx / r;
force[1] += factor * pdy / r;
force[2] += factor * pdz / r;
}
}
return force;
}
double r2 = pdx*pdx + pdy*pdy + pdz*pdz;
double r = std::sqrt(r2);
double factor = mass_*mass_ / (r2 + softening_*softening_);
Vector Node::multipole(Vector y, double ynorm) const {
Vector result = {0, 0, 0};
for (unsigned i=0; i<3; ++i) {
result[i] += totalMass_*mass_/(ynorm*ynorm) * y[i]/ynorm;
}
// calculate matrix Q
Matrix Q = {{0,0,0},{0,0,0},{0,0,0}};
for (unsigned i=0; i<3; ++i) {
for (unsigned j=0; j<3; ++j) {
// calculate the entry in Qij
double qij = 0;
// go through every local particle
for (unsigned c=0; c<8; ++c) {
for (auto k : childParticles_[c]) {
// calculate using the formula from page 17 of the pdf for lecture 7
qij += 3*(centerOfMass_[i] - positions_[3*k+i])*(centerOfMass_[j] - positions_[3*k+j]);
if (i!=j) continue;
for (unsigned l=0; l<3; ++l) {
qij -= (centerOfMass_[l] - positions_[3*k+l])*(centerOfMass_[l] - positions_[3*k+l]);
}
// add the force from each particle in the child
force[0] += factor * pdx / r;
force[1] += factor * pdy / r;
force[2] += factor * pdz / r;
}
}
}
}
} else {
// calculate force from this node (from the center of mass) with multipole expansion
// TODO: calculate Q and implement multipole expansion
double r2 = lambda*lambda;
double factor = totalMass_*mass_ / (r2 + softening_*softening_);
// y^T * Q * y
double yTQy = 0;
for (unsigned i=0; i<3; ++i) {
for (unsigned j=0; j<3; ++j) {
yTQy += y[i]*Q[i][j]*y[j];
}
}
// (Q + Q^T) * y
Vector QQTy = {0, 0, 0};
for (unsigned i=0; i<3; ++i) {
for (unsigned j=0; j<3; ++j) {
QQTy[i] += (Q[i][j] + Q[j][i])*y[j];
}
}
force[0] += factor * dx / lambda;
force[1] += factor * dy / lambda;
force[2] += factor * dz / lambda;
for (unsigned i=0; i<3; ++i) {
result[i] += 0.5*(-5 * yTQy / ynorm + QQTy[i] * y[i])/(ynorm * ynorm*ynorm*ynorm*ynorm);
}
return force;
return result;
}
void Node::updateMassAndCOM() {
......
lib/node.hpp
+ 3
− 1
View file @ f843e07d
......@@ -4,6 +4,7 @@
#include <vector>
using Vector = std::vector<double>;
using Matrix = std::vector<Vector>;
using List = std::vector<unsigned>;
class Node {
......@@ -24,6 +25,8 @@ public:
void update(List particles, unsigned nParticles);
// Calculate force on a particle recursively
Vector calculateForce(unsigned particle);
// Calculate the matrix Q
Vector multipole(Vector y, double ynorm) const;
// Update center of mass and total mass in the node
void updateMassAndCOM();
// Function that takes a particle index in the global particles_ array and returns the octant it is in relative to the center fo the current node
......@@ -48,7 +51,6 @@ private:
// Contains the indices of the particles that are in the children of this node for the particles_ array
List childParticles_[8];
unsigned nLocalParticles_;
double Q_[9];
double totalMass_;
const double theta0_ = 0.1;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment