Google Git
Sign in
gem5 / public / gem5-resources / 32b3c7eef2c389f18d9746d3ac5b50f72ff588d8 / . / src / halo-finder / src / BHForceTree.cxx
blob: b8c9c997064f99363f2470417d58c7b928b6110b [file] [log] [blame]
/*=========================================================================
Copyright (c) 2007, Los Alamos National Security, LLC
All rights reserved.
Copyright 2007. Los Alamos National Security, LLC.
This software was produced under U.S. Government contract DE-AC52-06NA25396
for Los Alamos National Laboratory (LANL), which is operated by
Los Alamos National Security, LLC for the U.S. Department of Energy.
The U.S. Government has rights to use, reproduce, and distribute this software.
NEITHER THE GOVERNMENT NOR LOS ALAMOS NATIONAL SECURITY, LLC MAKES ANY WARRANTY,
EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.
If software is modified to produce derivative works, such modified software
should be clearly marked, so as not to confuse it with the version available
from LANL.
Additionally, redistribution and use in source and binary forms, with or
without modification, are permitted provided that the following conditions
are met:
- Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation
and/or other materials provided with the distribution.
- Neither the name of Los Alamos National Security, LLC, Los Alamos National
Laboratory, LANL, the U.S. Government, nor the names of its contributors
may be used to endorse or promote products derived from this software
without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY LOS ALAMOS NATIONAL SECURITY, LLC AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL LOS ALAMOS NATIONAL SECURITY, LLC OR
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=========================================================================*/
#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <set>
#include <math.h>
#include <time.h>
#include "Timings.h"
#include "BHForceTree.h"
using namespace std;
/////////////////////////////////////////////////////////////////////////
//
// FParticle contains information about particles
//
/////////////////////////////////////////////////////////////////////////
FParticle::FParticle()
{
this->parent = -1;
this->nextNode = -1;
this->sibling = -1;
this->force = 0.0;
}
/////////////////////////////////////////////////////////////////////////
//
// FNode is a region of physical space divided into octants
//
/////////////////////////////////////////////////////////////////////////
FNode::FNode(POSVEL_T* minLoc, POSVEL_T* maxLoc)
{
for (int dim = 0; dim < DIMENSION; dim++) {
this->geoSide[dim] = maxLoc[dim] - minLoc[dim];
this->geoCenter[dim] = minLoc[dim] + this->geoSide[dim] * 0.5;
}
for (int i = 0; i < NUM_CHILDREN; i++)
this->u.child[i] = -1;
}
/////////////////////////////////////////////////////////////////////////
//
// FNode constructed from an octant of a parent node
//
/////////////////////////////////////////////////////////////////////////
FNode::FNode(FNode* parent, int oindx)
{
for (int dim = 0; dim < DIMENSION; dim++) {
this->geoSide[dim] = parent->geoSide[dim] * 0.5;
}
// Vary Z fastest when making octtree children
// If this changes must also change getChildIndex()
if (oindx & 4)
this->geoCenter[0] = parent->geoCenter[0] + this->geoSide[0] * 0.5;
else
this->geoCenter[0] = parent->geoCenter[0] - this->geoSide[0] * 0.5;
if (oindx & 2)
this->geoCenter[1] = parent->geoCenter[1] + this->geoSide[1] * 0.5;
else
this->geoCenter[1] = parent->geoCenter[1] - this->geoSide[1] * 0.5;
if (oindx & 1)
this->geoCenter[2] = parent->geoCenter[2] + this->geoSide[2] * 0.5;
else
this->geoCenter[2] = parent->geoCenter[2] - this->geoSide[2] * 0.5;
for (int i = 0; i < NUM_CHILDREN; i++)
this->u.child[i] = -1;
}
/////////////////////////////////////////////////////////////////////////
//
// Barnes Hut Tree
//
/////////////////////////////////////////////////////////////////////////
BHForceTree::BHForceTree(
POSVEL_T* minLoc,
POSVEL_T* maxLoc,
ID_T count,
POSVEL_T* xLoc,
POSVEL_T* yLoc,
POSVEL_T* zLoc,
POSVEL_T* xVel,
POSVEL_T* yVel,
POSVEL_T* zVel,
POSVEL_T* ms,
POSVEL_T avgMass)
{
// Extract the contiguous data block from a vector pointer
this->particleCount = count;
this->nodeOffset = this->particleCount;
this->xx = xLoc;
this->yy = yLoc;
this->zz = zLoc;
this->vx = xVel;
this->vy = yVel;
this->vz = zVel;
this->mass = ms;
this->particleMass = avgMass;
// Find the grid size of this chaining mesh
for (int dim = 0; dim < DIMENSION; dim++) {
this->minRange[dim] = minLoc[dim];
this->maxRange[dim] = maxLoc[dim];
}
this->boxSize = this->maxRange[0] - this->minRange[0];
//maybe change this to Newton's law or something
this->m_fl = new ForceLawNewton();
this->m_fcoeff = 1.0;
// Create the recursive BH tree from the particle locations
createBHForceTree();
// Thread the recursive tree turning it into an iterative tree
ID_T rootIndx = this->particleCount;
ID_T sibling = -1;
ID_T parent = -1;
ID_T lastIndx = -1;
POSVEL_T radius = 0.0;
threadBHForceTree(rootIndx, sibling, parent, &lastIndx, &radius);
}
BHForceTree::BHForceTree(
POSVEL_T* minLoc,
POSVEL_T* maxLoc,
ID_T count,
POSVEL_T* xLoc,
POSVEL_T* yLoc,
POSVEL_T* zLoc,
POSVEL_T* xVel,
POSVEL_T* yVel,
POSVEL_T* zVel,
POSVEL_T* ms,
POSVEL_T avgMass,
ForceLaw *fl,
float fcoeff)
{
// Extract the contiguous data block from a vector pointer
this->particleCount = count;
this->nodeOffset = this->particleCount;
this->xx = xLoc;
this->yy = yLoc;
this->zz = zLoc;
this->vx = xVel;
this->vy = yVel;
this->vz = zVel;
this->mass = ms;
this->particleMass = avgMass;
// Find the grid size of this chaining mesh
for (int dim = 0; dim < DIMENSION; dim++) {
this->minRange[dim] = minLoc[dim];
this->maxRange[dim] = maxLoc[dim];
}
this->boxSize = this->maxRange[0] - this->minRange[0];
this->m_fl = fl;
this->m_fcoeff = fcoeff;
// Create the recursive BH tree from the particle locations
createBHForceTree();
// Thread the recursive tree turning it into an iterative tree
ID_T rootIndx = this->particleCount;
ID_T sibling = -1;
ID_T parent = -1;
ID_T lastIndx = -1;
POSVEL_T radius = 0.0;
threadBHForceTree(rootIndx, sibling, parent, &lastIndx, &radius);
}
BHForceTree::~BHForceTree()
{
/* empty */
}
/////////////////////////////////////////////////////////////////////////
//
// Find the subhalos of the FOF halo using SUBFIND algorithm which
// requires subhalos to be locally overdense and self-bound
//
/////////////////////////////////////////////////////////////////////////
void BHForceTree::createBHForceTree()
{
// Create the FParticles
this->fParticle.resize(this->particleCount);
// Reserve a basic amount of space in the BH tree
this->fNode.reserve(this->particleCount/NUM_CHILDREN);
// Create the root node of the BH tree
FNode root(this->minRange, this->maxRange);
this->fNode.push_back(root);
ID_T nodeIndex = 0;
// Iterate on all particles placing them in the BH tree
// Child slots in the tree contain the index of the FParticle or
// the index of the FNode offset by the number of particles
// This is so we can use an integer instead of pointers to refer to objects
//
for (ID_T pindx = 0; pindx < this->particleCount; pindx++) {
// Start at root of tree for insertion of a new particle
// pindx is index into the halo particles where location is stored
// tindx is index into the BH tree nodes
// oindx is index into the octant of the tree node
ID_T tindx = 0;
int oindx = getChildIndex(&this->fNode[tindx], pindx);
while (this->fNode[tindx].u.child[oindx] != -1) {
// Child slot in tree contains another SPHNode so go there
if (this->fNode[tindx].u.child[oindx] > this->particleCount) {
tindx = this->fNode[tindx].u.child[oindx] - this->particleCount;
oindx = getChildIndex(&this->fNode[tindx], pindx);
}
// Otherwise there is a particle in the slot and we make a new FNode
else {
// Get the particle index of particle already in the node
ID_T pindx2 = this->fNode[tindx].u.child[oindx];
// First, check to make sure that this particle is not at the exact
// same location as the particle that is already there. If it is, then
// we'll double the mass of the existing particle and leave this one
// out.
if (this->xx[pindx2] == this->xx[pindx] &&
this->yy[pindx2] == this->yy[pindx] &&
this->zz[pindx2] == this->zz[pindx]) {
this->mass[pindx2] += this->mass[pindx];
goto next_particle;
}
// Make sure that the vector does not over allocate
if (this->fNode.capacity() == this->fNode.size()) {
this->fNode.reserve(this->fNode.capacity()
+ this->particleCount/NUM_CHILDREN);
}
FNode node(&this->fNode[tindx], oindx);
this->fNode.push_back(node);
nodeIndex++;
ID_T tindx2 = nodeIndex;
// Place the node that was sitting there already
int oindx2 = getChildIndex(&this->fNode[tindx2], pindx2);
this->fNode[tindx2].u.child[oindx2] = pindx2;
// Add the new SPHNode to the BHTree
this->fNode[tindx].u.child[oindx] = tindx2 + this->particleCount;
// Set to new node
tindx = tindx2;
oindx = getChildIndex(&this->fNode[tindx], pindx);
}
}
// Place the current particle in the BH tree
this->fNode[tindx].u.child[oindx] = pindx;
next_particle:;
}
this->nodeCount = this->fNode.size();
}
/////////////////////////////////////////////////////////////////////////
//
// Update the FNode vector by walking using a depth first recursion
// Set parent and sibling indices which can replace the child[8] already
// there, and supply extra information about center of mass, total particle
// mass and particle radius which is the distance from the center of mass
// to the furthest particle.
//
/////////////////////////////////////////////////////////////////////////
void BHForceTree::threadBHForceTree(
ID_T curIndx, // Current node/particle
ID_T sibling, // Sibling of current
ID_T parent, // Parent of current
ID_T* lastIndx, // Last node/particle
POSVEL_T* radius) // Needed to pass up partRadius
{
// Set the next index in the threading for node or particle
// Particles and nodes are threaded together so all are touched in iteration
if (*lastIndx >= 0) {
if (*lastIndx >= this->nodeOffset)
this->fNode[(*lastIndx - this->nodeOffset)].u.n.nextNode = curIndx;
else
this->fParticle[*lastIndx].nextNode = curIndx;
}
*lastIndx = curIndx;
// FParticle saves the parent and sibling FNode id
if (curIndx < this->nodeOffset) {
this->fParticle[curIndx].parent = parent;
this->fParticle[curIndx].sibling = sibling;
// FNode recurses on each of the children
} else {
ID_T child[NUM_CHILDREN];
FNode* curNode = &this->fNode[curIndx - this->nodeOffset];
// Store mass and center of mass for each child node or particle
POSVEL_T childMass[NUM_CHILDREN];
POSVEL_T childRadius[NUM_CHILDREN];
POSVEL_T childCenter[NUM_CHILDREN][DIMENSION];
for (int j = 0; j < NUM_CHILDREN; j++) {
child[j] = curNode->u.child[j];
childMass[j] = 0.0;
childRadius[j] = 0.0;
}
////////////////////////////////////////////////////////////////////////
//
// Recurse on each of the children, recording information on the way up
//
for (int j = 0; j < NUM_CHILDREN; j++) {
// Skip any children which contain neither a particle or node
ID_T childIndx, childIndxNext, nextSibling;
if ((childIndx = child[j]) >= 0) {
// Check for a sibling on this level or set to the next level up
int jj;
for (jj = j + 1; jj < NUM_CHILDREN; jj++)
if ((childIndxNext = child[jj]) >= 0)
break;
if (jj < NUM_CHILDREN)
nextSibling = childIndxNext;
else
nextSibling = -1;
// Recursion to child
// Since value of partRadius set in child is not necessarily the
// distance between center of mass and futhest child return it
threadBHForceTree(childIndx, nextSibling, curIndx, lastIndx, radius);
// Child is a node or a particle
if (childIndx >= this->nodeOffset) {
// FNode, gather mass and center of mass of all contained particles
FNode* childNode = &this->fNode[childIndx - this->nodeOffset];
childMass[j] = childNode->u.n.partMass;
childRadius[j] = *radius;
for (int dim = 0; dim < DIMENSION; dim++)
childCenter[j][dim] = childNode->u.n.partCenter[dim];
} else {
// FParticle, set mass and center of mass using particle location
childMass[j] = this->particleMass;
childRadius[j] = 0.0;
childCenter[j][0] = this->xx[childIndx];
childCenter[j][1] = this->yy[childIndx];
childCenter[j][2] = this->zz[childIndx];
}
}
}
////////////////////////////////////////////////////////////////////////
//
// Finished processing all children, collect information for this node
//
curNode->u.n.partMass = 0.0;
for (int dim = 0; dim < DIMENSION; dim++)
curNode->u.n.partCenter[dim] = 0.0;
// Collect total mass and total center of mass for all children
for (int j = 0; j < NUM_CHILDREN; j++) {
if (childMass[j] > 0) {
curNode->u.n.partMass += childMass[j];
for (int dim = 0; dim < DIMENSION; dim++)
curNode->u.n.partCenter[dim] +=
childCenter[j][dim] * childMass[j];
}
}
// Calculate center of mass for current node
if (curNode->u.n.partMass > 0.0) {
for (int dim = 0; dim < DIMENSION; dim++)
curNode->u.n.partCenter[dim] /= curNode->u.n.partMass;
} else {
for (int dim = 0; dim < DIMENSION; dim++)
curNode->u.n.partCenter[dim] = curNode->geoCenter[dim];
}
// First method for calculating particle radius
// Calculate the radius from node center of mass to furthest node corner
POSVEL_T partRadius1 = distanceToFarCorner(
curNode->u.n.partCenter[0],
curNode->u.n.partCenter[1],
curNode->u.n.partCenter[2],
curNode);
// Second method for calculating particle radius
// Calculate the radius from center of mass to furthest child
POSVEL_T partRadius2 = 0.0;
for (int j = 0; j < NUM_CHILDREN; j++) {
if (childMass[j] > 0.0) {
// Calculate the distance between this center of mass and that of child
POSVEL_T dist = distanceToCenterOfMass(childCenter[j][0],
childCenter[j][1],
childCenter[j][2],
curNode);
// Add in the particle radius of the child to get furthest point
dist += childRadius[j];
if (dist > partRadius2)
partRadius2 = dist;
}
}
// Used by parent of this node
*radius = partRadius2;
// Save the smaller of the two particle radii
curNode->u.n.partRadius = min(partRadius1, partRadius2);
// Set threaded structure for this node
curNode->u.n.sibling = sibling;
curNode->u.n.parent = parent;
}
}
/////////////////////////////////////////////////////////////////////////
//
// Particle i has location vector x_i (xx_i, yy_i, zz_i)
// Particle i has velocity vector v_i (vx_i, vy_i, vz_i)
//
// Node j has center of particle mass c_j (cx_j, cy_j, cz_j)
// Node j has total particle mass of M_j
// Node j has bounding radius of R_j
//
// Distance between particle i and node j is
// d_ji = fabs(c_j - x_i)
//
// Rule for updating
// v_i(t_1) = v_i(t_0) + alpha * Sum over j of f_ji(d_ji, Mj)
// where f_ji is the short range force over finite range r_f
// where alpha is some coeffient
// where Sum over j nodes is determined by a tree walk
//
// An opening angle is defined as
// theta_ji = (2 R_j) / d_ji
//
// This angle determines whether a node should be opened to a higher resolution
// or whether it can be used as is because it is small enough or far enough away
// This is determined by comparing to a passed in theta_0
//
// Three actions can occur for a node encountered on the walk
//
// 1. Node is too far away to contribute to force
// if d_ji - R_j > r_f
// or distance of x_i to nearest cornder of node > r_f
//
// 2. Node is close enough to contribute so check the opening angle
// if theta_ji > theta_0 follow nextNode to open this node to children
//
// 3. Node is close enough and theta_ji < theta_0
// calculate f_ji(d_ji, Mj) and update v_i
// follow the sibling link and not the nextNode link
//
// Force is calculated for each particle i by
// Starting at the root node and walking the entire tree collecting force
// Starting at the particle and walking up parents until a criteria is met
//
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
//
// Short range force full N^2 calculation
//
/////////////////////////////////////////////////////////////////////////
void BHForceTree::treeForceN2(
POSVEL_T critRadius) // Radius of furthest point
{
POSVEL_T critRadius2 = critRadius * critRadius;
POTENTIAL_T* force = new POTENTIAL_T[this->particleCount];
for (int i = 0; i < this->particleCount; i++)
force[i] = 0.0;
// First particle in halo to calculate force on
for (int p = 0; p < this->particleCount; p++) {
// Next particle in halo force loop
for (int q = p+1; q < this->particleCount; q++) {
POSVEL_T dx = (POSVEL_T) fabs(this->xx[p] - this->xx[q]);
POSVEL_T dy = (POSVEL_T) fabs(this->yy[p] - this->yy[q]);
POSVEL_T dz = (POSVEL_T) fabs(this->zz[p] - this->zz[q]);
POSVEL_T r2 = dx * dx + dy * dy + dz * dz;
if (r2 != 0.0 && r2 < critRadius2) {
force[p] -= (this->mass[q] / r2);
force[q] -= (this->mass[p] / r2);
}
}
}
for (int p = 0; p < this->particleCount; p++) {
this->fParticle[p].force = force[p];
}
delete [] force;
}
/////////////////////////////////////////////////////////////////////////
//
// Short range gravity calculation for group of particles in a node
// Walk down the tree from the root until reaching node with less than the
// maximum number of particles in a group. Create an interaction list that
// will work for all particles and calculate force. For particles within
// the group the calculation will be n^2. For nodes outside the group
// decisions are made on whether to include the node or ignore it, or
// to accept it or open it.
//
/////////////////////////////////////////////////////////////////////////
void BHForceTree::treeForceGroup(
POSVEL_T bhAngle, // Open node to examine children
POSVEL_T critRadius, // Accept or ignore node not opened
int minGroup, // Minimum particles in a group
int maxGroup) // Maximum particles in a group
{
ID_T root = this->particleCount;
POSVEL_T maxMass = maxGroup * this->particleMass;
POSVEL_T minMass = minGroup * this->particleMass;
walkTreeGroup(root, minMass, maxMass, bhAngle, critRadius);
}
/////////////////////////////////////////////////////////////////////////
//
// Walk the tree in search of nodes which are less than the maximum
// number of particles to constitute a group. All particles in the group
// will be treated together with the n^2 force calculated between members
// of the group and then having an interaction list applied. The group
// may consist of other nodes and particles and so the recursive descent
// will be needed to find all particles in the group.
//
/////////////////////////////////////////////////////////////////////////
void BHForceTree::walkTreeGroup(
ID_T curId, // Particle to calculate force on
POSVEL_T minMass, // Minimum mass for a group
POSVEL_T maxMass, // Maximum mass for a group
POSVEL_T bhAngle, // Open node to examine children
POSVEL_T critRadius) // Accept or ignore node not opened
{
if (curId < this->nodeOffset) {
// Current object is a particle
vector<ID_T>* partInteract = new vector<ID_T>;
vector<ID_T>* nodeInteract = new vector<ID_T>;
createParticleInteractList(curId, bhAngle, critRadius,
partInteract, nodeInteract);
this->fParticle[curId].force = forceCalculation(
curId, partInteract, nodeInteract);
delete partInteract;
delete nodeInteract;
}
else {
// Current object is a node
ID_T child = this->fNode[curId - this->nodeOffset].u.n.nextNode;
while (child != -1) {
if (child < this->nodeOffset) {
// Child is a particle
vector<ID_T>* partInteract = new vector<ID_T>;
vector<ID_T>* nodeInteract = new vector<ID_T>;
createParticleInteractList(child, bhAngle, critRadius,
partInteract, nodeInteract);
this->fParticle[child].force = forceCalculation(
child, partInteract, nodeInteract);
child = this->fParticle[child].sibling;
delete partInteract;
delete nodeInteract;
}
else {
// Child is a node
FNode* childNode = &this->fNode[child - this->nodeOffset];
if (childNode->u.n.partMass < maxMass &&
childNode->u.n.partRadius < (critRadius * 0.5)) {
// If the group is too small it can't function as a group
// so run the topdown method on those particles
if (childNode->u.n.partMass < minMass) {
// Collect particles in subgroup
vector<ID_T>* particles = new vector<ID_T>;
collectParticles(child, particles);
int count = particles->size();
for (int i = 0; i < count; i++) {
ID_T pId = (*particles)[i];
treeForceGadgetTopDown(pId, bhAngle, critRadius);
}
}
else {
vector<ID_T>* partInteract = new vector<ID_T>;
vector<ID_T>* nodeInteract = new vector<ID_T>;
createNodeInteractList(child, bhAngle, critRadius,
partInteract, nodeInteract);
forceCalculationGroup(child, bhAngle, critRadius,
partInteract, nodeInteract);
delete partInteract;
delete nodeInteract;
}
}
else {
walkTreeGroup(child, minMass, maxMass, bhAngle, critRadius);
}
child = this->fNode[child - this->nodeOffset].u.n.sibling;
}
}
}
}
/////////////////////////////////////////////////////////////////////////
//
// Create the interaction list for the particle starting at root
//
/////////////////////////////////////////////////////////////////////////
void BHForceTree::createParticleInteractList(
ID_T p, // Particle to calculate force on
POSVEL_T bhAngle, // Open node to examine children
POSVEL_T critRadius , // Accept or ignore node not opened
vector<ID_T>* partInteract,
vector<ID_T>* nodeInteract)
{
POSVEL_T dx, dy, dz, r, partRadius, distToNearPoint;
POSVEL_T pos_x = this->xx[p];
POSVEL_T pos_y = this->yy[p];
POSVEL_T pos_z = this->zz[p];
// Follow thread through tree from root choosing nodes and particles
// which will contribute to the force of the given particle
ID_T root = this->particleCount;
ID_T index = root;
while (index >= 0) {
if (index < this->nodeOffset) {
// Particle
dx = this->xx[index] - pos_x;
dy = this->yy[index] - pos_y;
dz = this->zz[index] - pos_z;
r = sqrt(dx * dx + dy * dy + dz * dz);
if (r < critRadius) {
partInteract->push_back(index);
}
index = this->fParticle[index].nextNode;
}
else {
// Node
FNode* curNode = &this->fNode[index - this->nodeOffset];
partRadius = curNode->u.n.partRadius;
distToNearPoint = distanceToNearestPoint(pos_x, pos_y, pos_z, curNode);
dx = curNode->u.n.partCenter[0] - pos_x;
dy = curNode->u.n.partCenter[1] - pos_y;
dz = curNode->u.n.partCenter[2] - pos_z;
r = sqrt(dx * dx + dy * dy + dz * dz);
// Node is ignored if it is too far away from the particle
// Distance from particle to particle radius exceeds critical radius
// Distance from particle to nearest side of node exceeds critical radius
if ((r - partRadius) > critRadius || distToNearPoint > critRadius) {
// Ignore node, move on to sibling of this node
index = curNode->u.n.sibling;
// If there is no sibling go up a level until we find a node
ID_T parent = curNode->u.n.parent;
while (index == -1 && parent != -1 && parent != root) {
index = this->fNode[parent - this->nodeOffset].u.n.sibling;
parent = this->fNode[parent - this->nodeOffset].u.n.parent;
}
}
else {
if (2*partRadius > (r * bhAngle)) {
// Open node
index = curNode->u.n.nextNode;
} else {
// Accept
nodeInteract->push_back(index);
index = curNode->u.n.sibling;
// If there is no sibling go up a level until we find a node
ID_T parent = curNode->u.n.parent;
while (index == -1 && parent != -1 && parent != root) {
index = this->fNode[parent - this->nodeOffset].u.n.sibling;
parent = this->fNode[parent - this->nodeOffset].u.n.parent;
}
}
}
}
}
}
/////////////////////////////////////////////////////////////////////////
//
// Create the interaction list for the node starting at root
// must test for acceptance based on a radius from center of mass to
// furthest particle to make sure it is most inclusive.
// Make sure my definition of partRadius does this
//
/////////////////////////////////////////////////////////////////////////
void BHForceTree::createNodeInteractList(
ID_T node,
POSVEL_T bhAngle, // Open node to examine children
POSVEL_T critRadius, // Accept or ignore node not opened
vector<ID_T>* partInteract,
vector<ID_T>* nodeInteract)
{
POSVEL_T dx, dy, dz, r, partRadius, distToNearPoint;
FNode* curNode = &this->fNode[node - this->nodeOffset];
POSVEL_T pos_x = curNode->u.n.partCenter[0];
POSVEL_T pos_y = curNode->u.n.partCenter[1];
POSVEL_T pos_z = curNode->u.n.partCenter[2];
// Follow thread through tree from root choosing nodes and particles
// which will contribute to the force of the given particle
ID_T root = this->particleCount;
ID_T index = root;
while (index >= 0) {
if (index < this->nodeOffset) {
// Particle
dx = this->xx[index] - pos_x;
dy = this->yy[index] - pos_y;
dz = this->zz[index] - pos_z;
r = sqrt(dx * dx + dy * dy + dz * dz);
if (r < critRadius) {
partInteract->push_back(index);
}
index = this->fParticle[index].nextNode;
}
else {
// Node
FNode* childNode = &this->fNode[index - this->nodeOffset];
// If the child is the node we are building the list for skip
if (childNode != curNode) {
partRadius = childNode->u.n.partRadius;
distToNearPoint = distanceToNearestPoint(
pos_x, pos_y, pos_z, childNode);
dx = childNode->u.n.partCenter[0] - pos_x;
dy = childNode->u.n.partCenter[1] - pos_y;
dz = childNode->u.n.partCenter[2] - pos_z;
r = sqrt(dx * dx + dy * dy + dz * dz);
if ((r - partRadius) > critRadius || distToNearPoint > critRadius) {
// Ignore node, move on to sibling of this node
index = childNode->u.n.sibling;
// If there is no sibling go up a level until we find a node
ID_T parent = childNode->u.n.parent;
while (index == -1 && parent != -1 && parent != root) {
index = this->fNode[parent - this->nodeOffset].u.n.sibling;
parent = this->fNode[parent - this->nodeOffset].u.n.parent;
}
}
else {
if (2*partRadius > (r * bhAngle)) {
// Open node
index = childNode->u.n.nextNode;
} else {
// Accept
nodeInteract->push_back(index);
index = childNode->u.n.sibling;
// If there is no sibling go up a level until we find a node
ID_T parent = childNode->u.n.parent;
while (index == -1 && parent != -1 && parent != root) {
index = this->fNode[parent - this->nodeOffset].u.n.sibling;
parent = this->fNode[parent - this->nodeOffset].u.n.parent;
}
}
}
}
else {
index = childNode->u.n.sibling;
ID_T parent = childNode->u.n.parent;
while (index == -1 && parent != -1 && parent != root) {
index = this->fNode[parent-nodeOffset].u.n.sibling;
parent = this->fNode[parent-nodeOffset].u.n.parent;
}
}
}
}
}
/////////////////////////////////////////////////////////////////////////
//
// Force calculation for a group of particles
// Force is calculated between every pair of particles in the group
// Interaction lists are applied to every particle in the group
// Tree walk will hav to continue from this node to locate all the
// particles which might be in subnodes of this node.
//
/////////////////////////////////////////////////////////////////////////
void BHForceTree::forceCalculationGroup(
ID_T node,
POSVEL_T bhAngle, // Open or accepc
POSVEL_T critRadius, // Accept or ignore node not opened
vector<ID_T>* partInteract,
vector<ID_T>* nodeInteract)
{
// Collect all particles in the tree from this node downwards
vector<ID_T>* particles = new vector<ID_T>;
collectParticles(node, particles);
int count = particles->size();
// Process each particle against all others in the group
// Use the minimumPotential() code to make this n^2/2 for upper triangular
// Arrange in an upper triangular grid to save computation
POTENTIAL_T* force = new POTENTIAL_T[count];
for (int i = 0; i < count; i++)
force[i] = 0.0;
// First particle in halo to calculate force on
for (int p = 0; p < count; p++) {
// Next particle in halo force loop
for (int q = p+1; q < count; q++) {
ID_T pId = (*particles)[p];
ID_T qId = (*particles)[q];
POSVEL_T dx = (POSVEL_T) fabs(this->xx[pId] - this->xx[qId]);
POSVEL_T dy = (POSVEL_T) fabs(this->yy[pId] - this->yy[qId]);
POSVEL_T dz = (POSVEL_T) fabs(this->zz[pId] - this->zz[qId]);
POSVEL_T r2 = dx * dx + dy * dy + dz * dz;
if (r2 != 0.0) {
force[p] -= (this->mass[qId] / r2);
force[q] -= (this->mass[pId] / r2);
}
}
}
// Process each particle against the interaction lists
// Node interact list was created using the node this particle is in
// so it may need to be adjusted first
for (int p = 0; p < count; p++) {
ID_T pId = (*particles)[p];
POSVEL_T value =
forceCalculationParticle(pId, critRadius,
partInteract, nodeInteract);
force[p] += value;
this->fParticle[pId].force = force[p];
}
delete particles;
delete [] force;
}
/////////////////////////////////////////////////////////////////////////
//
// Short range force calculation
// Potential is calculated and is used to determine the acceleration of
// the particle. Acceleration is applied to the current velocity to
// produce the velocity at the next time step.
//
/////////////////////////////////////////////////////////////////////////
POSVEL_T BHForceTree::forceCalculationParticle(
ID_T p0, // Target particle index
POSVEL_T critRadius,
vector<ID_T>* partInteract, // Particles acting on p
vector<ID_T>* nodeInteract) // Nodes acting on p
{
POSVEL_T accel[DIMENSION];
POSVEL_T phi = 0.0;
POSVEL_T pos0_x = this->xx[p0];
POSVEL_T pos0_y = this->yy[p0];
POSVEL_T pos0_z = this->zz[p0];
for (int dim = 0; dim < DIMENSION; dim++)
accel[dim] = 0.0;
int numberOfNodes = (int) nodeInteract->size();
int numberOfParticles = (int) partInteract->size();
// Particles contributing to the force use location and mass of one particle
for (int p = 0; p < numberOfParticles; p++) {
ID_T particle = (*partInteract)[p];
if (p0 != particle) {
POSVEL_T dx = this->xx[particle] - pos0_x;
POSVEL_T dy = this->yy[particle] - pos0_y;
POSVEL_T dz = this->zz[particle] - pos0_z;
POSVEL_T r2 = dx * dx + dy * dy + dz * dz;
POSVEL_T r = sqrt(r2);
if (r < critRadius) {
POSVEL_T f_over_r = this->mass[particle] * m_fl->f_over_r(r2);
//POSVEL_T f_over_r = this->mass[particle] / r2;
phi -= f_over_r;
accel[0] += dx * f_over_r * m_fcoeff;
accel[1] += dy * f_over_r * m_fcoeff;
accel[2] += dz * f_over_r * m_fcoeff;
this->vx[p0] += dx * f_over_r * m_fcoeff;
this->vy[p0] += dy * f_over_r * m_fcoeff;
this->vz[p0] += dz * f_over_r * m_fcoeff;
}
}
}
// Nodes contributing to force use center of mass and total particle mass
for (int n = 0; n < numberOfNodes; n++) {
FNode* node = &this->fNode[(*nodeInteract)[n] - this->nodeOffset];
POSVEL_T partRadius = node->u.n.partRadius;
POSVEL_T distToNearPoint = distanceToNearestPoint(
pos0_x, pos0_y, pos0_z, node);
POSVEL_T dx = node->u.n.partCenter[0] - pos0_x;
POSVEL_T dy = node->u.n.partCenter[1] - pos0_y;
POSVEL_T dz = node->u.n.partCenter[2] - pos0_z;
POSVEL_T r2 = dx * dx + dy * dy + dz * dz;
POSVEL_T r = sqrt(r2);
if ((r - partRadius) > critRadius || distToNearPoint > critRadius) {
// Ignore
} else {
POSVEL_T f_over_r = node->u.n.partMass * m_fl->f_over_r(r2);
//POSVEL_T f_over_r = node->u.n.partMass / r2;
phi -= f_over_r;
accel[0] += dx * f_over_r * m_fcoeff;
accel[1] += dy * f_over_r * m_fcoeff;
accel[2] += dz * f_over_r * m_fcoeff;
this->vx[p0] += dx * f_over_r * m_fcoeff;
this->vy[p0] += dy * f_over_r * m_fcoeff;
this->vz[p0] += dz * f_over_r * m_fcoeff;
}
}
return phi;
}
/////////////////////////////////////////////////////////////////////////
//
// Collect all particle ids from this node downwards in tree
//
/////////////////////////////////////////////////////////////////////////
void BHForceTree::collectParticles(ID_T curId, vector<ID_T>* particles)
{
FNode* curNode = &this->fNode[curId - this->nodeOffset];
ID_T child = curNode->u.n.nextNode;
while (child != -1) {
if (child < this->nodeOffset) {
particles->push_back(child);
child = this->fParticle[child].sibling;
} else {
collectParticles(child, particles);
child = this->fNode[child - this->nodeOffset].u.n.sibling;
}
}
}
/////////////////////////////////////////////////////////////////////////
//
// Short range gravity calculation for a single particle
// Starting with the root and following threads and siblings makes decisions
// about which nodes are opened, accepted or ignored
//
/////////////////////////////////////////////////////////////////////////
void BHForceTree::treeForceGadgetTopDown(
ID_T p, // Particle to calculate force on
POSVEL_T bhAngle, // Open node to examine children
POSVEL_T critRadius) // Accept or ignore node not opened
{
// Keep vectors for now but eventually force can be accumulated
// on all particles and on nodes that have been accepted
vector<ID_T>* partInteract = new vector<ID_T>;
vector<ID_T>* nodeInteract = new vector<ID_T>;
POSVEL_T dx, dy, dz, r, partRadius, distToNearPoint;
POSVEL_T pos_x = this->xx[p];
POSVEL_T pos_y = this->yy[p];
POSVEL_T pos_z = this->zz[p];
// Follow thread through tree from root choosing nodes and particles
// which will contribute to the force of the given particle
ID_T root = this->particleCount;
ID_T index = root;
while (index >= 0) {
if (index < this->nodeOffset) {
// Particle
dx = this->xx[index] - pos_x;
dy = this->yy[index] - pos_y;
dz = this->zz[index] - pos_z;
r = sqrt(dx * dx + dy * dy + dz * dz);
if (r < critRadius) {
partInteract->push_back(index);
}
index = this->fParticle[index].nextNode;
} else {
// Node
FNode* curNode = &this->fNode[index - this->nodeOffset];
partRadius = curNode->u.n.partRadius;
distToNearPoint = distanceToNearestPoint(pos_x, pos_y, pos_z, curNode);
dx = curNode->u.n.partCenter[0] - pos_x;
dy = curNode->u.n.partCenter[1] - pos_y;
dz = curNode->u.n.partCenter[2] - pos_z;
r = sqrt(dx * dx + dy * dy + dz * dz);
// Node is ignored if it is too far away from the particle
// Distance from particle to particle radius exceeds critical radius
// Distance from particle to nearest side of node exceeds critical radius
if ((r - partRadius) > critRadius || distToNearPoint > critRadius) {
// Ignore node, move on to sibling of this node
index = curNode->u.n.sibling;
// If there is no sibling go up a level until we find a node
ID_T parent = curNode->u.n.parent;
while (index == -1 && parent != -1 && parent != root) {
index = this->fNode[parent - this->nodeOffset].u.n.sibling;
parent = this->fNode[parent - this->nodeOffset].u.n.parent;
}
} else {
if (2*partRadius > (r * bhAngle)) {
// Open node, move on to first child
index = curNode->u.n.nextNode;
} else {
// Accept node, add to interact list, move on to sibling
nodeInteract->push_back(index);
index = curNode->u.n.sibling;
// If there is no sibling go up a level until we find a node
ID_T parent = curNode->u.n.parent;
while (index == -1 && parent != -1 && parent != root) {
index = this->fNode[parent-nodeOffset].u.n.sibling;
parent = this->fNode[parent-nodeOffset].u.n.parent;
}
}
}
}
}
// Force calculation for this particle
this->fParticle[p].force =
forceCalculation(p, partInteract, nodeInteract);
delete partInteract;
delete nodeInteract;
}
void BHForceTree::treeForceGadgetTopDownFast(
ID_T p, // Particle to calculate force on
POSVEL_T bhAngle, // Open node to examine children
POSVEL_T critRadius) // Accept or ignore node not opened
{
vector<POSVEL_T>* xInteract = new vector<POSVEL_T>;
vector<POSVEL_T>* yInteract = new vector<POSVEL_T>;
vector<POSVEL_T>* zInteract = new vector<POSVEL_T>;
vector<POSVEL_T>* mInteract = new vector<POSVEL_T>;
POSVEL_T dx, dy, dz, r, partRadius, distToNearPoint;
POSVEL_T pos_x = this->xx[p];
POSVEL_T pos_y = this->yy[p];
POSVEL_T pos_z = this->zz[p];
// Follow thread through tree from root choosing nodes and particles
// which will contribute to the force of the given particle
ID_T root = this->particleCount;
ID_T index = root;
while (index >= 0) {
if (index < this->nodeOffset) {
// Particle
dx = this->xx[index] - pos_x;
dy = this->yy[index] - pos_y;
dz = this->zz[index] - pos_z;
r = sqrt(dx * dx + dy * dy + dz * dz);
if (r < critRadius && p != index) {
xInteract->push_back(this->xx[index]);
yInteract->push_back(this->yy[index]);
zInteract->push_back(this->zz[index]);
mInteract->push_back(this->mass[index]);
}
index = this->fParticle[index].nextNode;
} else {
// Node
FNode* curNode = &this->fNode[index - this->nodeOffset];
partRadius = curNode->u.n.partRadius;
distToNearPoint = distanceToNearestPoint(pos_x, pos_y, pos_z, curNode);
dx = curNode->u.n.partCenter[0] - pos_x;
dy = curNode->u.n.partCenter[1] - pos_y;
dz = curNode->u.n.partCenter[2] - pos_z;
r = sqrt(dx * dx + dy * dy + dz * dz);
// Node is ignored if it is too far away from the particle
// Distance from particle to particle radius exceeds critical radius
// Distance from particle to nearest side of node exceeds critical radius
if ((r - partRadius) > critRadius || distToNearPoint > critRadius) {
// Ignore node, move on to sibling of this node
index = curNode->u.n.sibling;
// If there is no sibling go up a level until we find a node
ID_T parent = curNode->u.n.parent;
while (index == -1 && parent != -1 && parent != root) {
index = this->fNode[parent - this->nodeOffset].u.n.sibling;
parent = this->fNode[parent - this->nodeOffset].u.n.parent;
}
} else {
if (2*partRadius > (r * bhAngle)) {
// Open node, move on to first child
index = curNode->u.n.nextNode;
} else {
// Accept node, add to interact list, move on to sibling
xInteract->push_back(curNode->u.n.partCenter[0]);
yInteract->push_back(curNode->u.n.partCenter[1]);
zInteract->push_back(curNode->u.n.partCenter[2]);
mInteract->push_back(curNode->u.n.partMass);
index = curNode->u.n.sibling;
// If there is no sibling go up a level until we find a node
ID_T parent = curNode->u.n.parent;
while (index == -1 && parent != -1 && parent != root) {
index = this->fNode[parent-nodeOffset].u.n.sibling;
parent = this->fNode[parent-nodeOffset].u.n.parent;
}
}
}
}
}
// Force calculation for this particle
this->fParticle[p].force =
forceCalculationFast(p, xInteract, yInteract, zInteract, mInteract);
delete xInteract;
delete yInteract;
delete zInteract;
delete mInteract;
}
void BHForceTree::treeForceGadgetTopDownFast2(
ID_T p, // Particle to calculate force on
POSVEL_T bhAngle, // Open node to examine children
POSVEL_T critRadius, // Accept or ignore node not opened
vector<POSVEL_T>* xInteract,
vector<POSVEL_T>* yInteract,
vector<POSVEL_T>* zInteract,
vector<POSVEL_T>* mInteract,
double *timeWalk,
double *timeEval)
{
POSVEL_T dx, dy, dz, r, partRadius, distToNearPoint;
POSVEL_T pos_x = this->xx[p];
POSVEL_T pos_y = this->yy[p];
POSVEL_T pos_z = this->zz[p];
// Follow thread through tree from root choosing nodes and particles
// which will contribute to the force of the given particle
ID_T root = this->particleCount;
ID_T index = root;
clock_t start, end;
start = clock();
while (index >= 0) {
if (index < this->nodeOffset) {
// Particle
dx = this->xx[index] - pos_x;
dy = this->yy[index] - pos_y;
dz = this->zz[index] - pos_z;
r = sqrt(dx * dx + dy * dy + dz * dz);
if (r < critRadius && p != index) {
xInteract->push_back(this->xx[index]);
yInteract->push_back(this->yy[index]);
zInteract->push_back(this->zz[index]);
mInteract->push_back(this->mass[index]);
}
index = this->fParticle[index].nextNode;
} else {
// Node
FNode* curNode = &this->fNode[index - this->nodeOffset];
partRadius = curNode->u.n.partRadius;
distToNearPoint = distanceToNearestPoint(pos_x, pos_y, pos_z, curNode);
dx = curNode->u.n.partCenter[0] - pos_x;
dy = curNode->u.n.partCenter[1] - pos_y;
dz = curNode->u.n.partCenter[2] - pos_z;
r = sqrt(dx * dx + dy * dy + dz * dz);
// Node is ignored if it is too far away from the particle
// Distance from particle to particle radius exceeds critical radius
// Distance from particle to nearest side of node exceeds critical radius
if ((r - partRadius) > critRadius || distToNearPoint > critRadius) {
// Ignore node, move on to sibling of this node
index = curNode->u.n.sibling;
// If there is no sibling go up a level until we find a node
ID_T parent = curNode->u.n.parent;
while (index == -1 && parent != -1 && parent != root) {
index = this->fNode[parent - this->nodeOffset].u.n.sibling;
parent = this->fNode[parent - this->nodeOffset].u.n.parent;
}
} else {
if (2*partRadius > (r * bhAngle)) {
// Open node, move on to first child
index = curNode->u.n.nextNode;
} else {
// Accept node, add to interact list, move on to sibling
xInteract->push_back(curNode->u.n.partCenter[0]);
yInteract->push_back(curNode->u.n.partCenter[1]);
zInteract->push_back(curNode->u.n.partCenter[2]);
mInteract->push_back(curNode->u.n.partMass);
index = curNode->u.n.sibling;
// If there is no sibling go up a level until we find a node
ID_T parent = curNode->u.n.parent;
while (index == -1 && parent != -1 && parent != root) {
index = this->fNode[parent-nodeOffset].u.n.sibling;
parent = this->fNode[parent-nodeOffset].u.n.parent;
}
}
}
}
}
end = clock();
*timeWalk = 1.0*(end-start)/CLOCKS_PER_SEC;
// Force calculation for this particle
start = clock();
this->fParticle[p].force =
forceCalculationFast(p, xInteract, yInteract, zInteract, mInteract);
end = clock();
*timeEval = 1.0*(end-start)/CLOCKS_PER_SEC;
xInteract->clear();
yInteract->clear();
zInteract->clear();
mInteract->clear();
}
/////////////////////////////////////////////////////////////////////////
//
// Short range gravity calculation for a single particle
// Starting with the particle walk up the parents processing siblings
// by testing particles and by opening, accepting or ignoring nodes.
// Stop moving up parents when the nearest side is beyond critical radius.
//
/////////////////////////////////////////////////////////////////////////
void BHForceTree::treeForceGadgetBottomUp(
ID_T p, // Particle to calculate force on
POSVEL_T bhAngle, // Open node to examine children
POSVEL_T critRadius) // Accept or ignore node not opened
{
// Collect into interaction lists
vector<ID_T>* partInteract = new vector<ID_T>;
vector<ID_T>* nodeInteract = new vector<ID_T>;
// Location of particle
POSVEL_T dx, dy, dz, r, partRadius;
POSVEL_T pos_x = this->xx[p];
POSVEL_T pos_y = this->yy[p];
POSVEL_T pos_z = this->zz[p];
ID_T curId = p;
ID_T parent = this->fParticle[curId].parent;
while (parent != -1) {
ID_T child = this->fNode[parent - this->nodeOffset].u.n.nextNode;
while (child != -1) {
if (child != curId) {
if (child < this->nodeOffset) {
// Particle
dx = this->xx[child] - pos_x;
dy = this->yy[child] - pos_y;
dz = this->zz[child] - pos_z;
r = sqrt(dx * dx + dy * dy + dz * dz);
if (r < critRadius) {
partInteract->push_back(child);
}
child = this->fParticle[child].sibling;
}
else {
// Node
FNode* childNode = &this->fNode[child - this->nodeOffset];
partRadius = childNode->u.n.partRadius;
POSVEL_T distToNearPoint =
distanceToNearestPoint(pos_x, pos_y, pos_z, childNode);
dx = childNode->u.n.partCenter[0] - pos_x;
dy = childNode->u.n.partCenter[1] - pos_y;
dz = childNode->u.n.partCenter[2] - pos_z;
r = sqrt(dx * dx + dy * dy + dz * dz);
// Check for ignore of node first
if ((r - partRadius) > critRadius || distToNearPoint > critRadius) {
// Ignore
} else
if (2*partRadius < (r * bhAngle)) {
// Accept
nodeInteract->push_back(child);
} else {
// Open node
recurseOpenNode(childNode, pos_x, pos_y, pos_z,
bhAngle, critRadius,
partInteract, nodeInteract);
}
child = this->fNode[child - this->nodeOffset].u.n.sibling;
}
}
else {
if (curId < this->nodeOffset)
child = this->fParticle[curId].sibling;
else
child = this->fNode[child - this->nodeOffset].u.n.sibling;
}
}
curId = parent;
parent = this->fNode[parent - this->nodeOffset].u.n.parent;
}
// Force calculation for this particle
this->fParticle[p].force =
forceCalculation(p, partInteract, nodeInteract);
delete partInteract;
delete nodeInteract;
}
/////////////////////////////////////////////////////////////////////////
//
// Open this node recursively adding accepted nodes and particles
// to the interact list
//
/////////////////////////////////////////////////////////////////////////
void BHForceTree::recurseOpenNode(
FNode* curNode,
POSVEL_T pos_x,
POSVEL_T pos_y,
POSVEL_T pos_z,
POSVEL_T bhAngle, // Open node to examine children
POSVEL_T critRadius, // Accept or ignore node not opened
vector<ID_T>* partInteract,
vector<ID_T>* nodeInteract)
{
POSVEL_T dx, dy, dz, r, partRadius;
ID_T child = curNode->u.n.nextNode;
while (child != -1) {
if (child < this->nodeOffset) {
// Particle
dx = this->xx[child] - pos_x;
dy = this->yy[child] - pos_y;
dz = this->zz[child] - pos_z;
r = sqrt(dx * dx + dy * dy + dz * dz);
if (r < critRadius) {
partInteract->push_back(child);
}
child = this->fParticle[child].sibling;
} else {
FNode* childNode = &this->fNode[child - this->nodeOffset];
partRadius = childNode->u.n.partRadius;
POSVEL_T distToNearPoint =
distanceToNearestPoint(pos_x, pos_y, pos_z, childNode);
dx = childNode->u.n.partCenter[0] - pos_x;
dy = childNode->u.n.partCenter[1] - pos_y;
dz = childNode->u.n.partCenter[2] - pos_z;
r = sqrt(dx * dx + dy * dy + dz * dz);
// Check for ignore of node first
if ((r - partRadius) > critRadius || distToNearPoint > critRadius) {
// Ignore
} else
if (2*partRadius < (r * bhAngle)) {
// Accept
nodeInteract->push_back(child);
} else {
// Open node
recurseOpenNode(childNode, pos_x, pos_y, pos_z,
bhAngle, critRadius,
partInteract, nodeInteract);
}
child = this->fNode[child - this->nodeOffset].u.n.sibling;
}
}
}
/////////////////////////////////////////////////////////////////////////
//
// Short range gravity calculation for every particle in the tree
// Recurses through the tree saving previous work for reuse when popping
// out of recursion. Based on Barnes treecode.
//
/////////////////////////////////////////////////////////////////////////
void BHForceTree::treeForceBarnesAdjust(
POSVEL_T bhAngle, // Open node to examine children
POSVEL_T critRadius) // Accept or ignore node not opened
{
ID_T root = this->particleCount;
vector<ID_T>* active = new vector<ID_T>;
vector<ID_T>* partInteract = new vector<ID_T>;
vector<ID_T>* nodeInteract = new vector<ID_T>;
active->push_back(root);
// Walk uses opening angle, critical radius for open, accept and ignore
walkTreeBarnesAdjust(active, partInteract, nodeInteract,
root, bhAngle, critRadius);
delete active;
delete partInteract;
delete nodeInteract;
}
///////////////////////////////////////////////////////////////////////////
//
// Walk the BH tree for the given particle or node (identifier curId)
// Recursion starts with a new active list which will contain particles
// and nodes which possibly will contribute to the force on a particle.
// Particles on the active list will always be chosen for the interact list.
// Nodes on the active list may be OPENED if they are close enough
// or ACCEPTED and used in summary if they are within a critical radius
// and IGNORED otherwise. Nodes that are opened
// have all their children (particles or nodes) added to the active list.
//
// After the children are added a new level of recursion starts by
// calculating a new size for that level, starting a fresh active list
// and building on the current interact lists.
//
// Recursion continues until the active list has been completely processed.
// When a level of recursion is complete the active list is destroyed
// and new items put on the interact lists are popped off.
//
// The advantage to this method is that items in the interaction list may
// not need to be processed again when we are doing the low levels of the tree.
//
/////////////////////////////////////////////////////////////////////////
void BHForceTree::walkTreeBarnesAdjust(
vector<ID_T>* curActive, // nodes to be acted on
vector<ID_T>* partInteract, // particles for force
vector<ID_T>* nodeInteract, // nodes for force
ID_T curId, // current particle or node
POSVEL_T bhAngle, // open node
POSVEL_T critRadius) // accept or ignore node
{
POSVEL_T dx, dy, dz, r, partRadius, distToNearPoint;
// Current active list
int begIndx = 0;
int endIndx = curActive->size();
// Construct active list for each recursion
vector<ID_T>* newActive = new vector<ID_T>;
// Set the location for the particle or node for the walk
POSVEL_T pos_x, pos_y, pos_z;
if (curId < this->nodeOffset) {
pos_x = this->xx[curId];
pos_y = this->yy[curId];
pos_z = this->zz[curId];
} else {
FNode* curNode = &this->fNode[curId - this->nodeOffset];
pos_x = curNode->u.n.partCenter[0];
pos_y = curNode->u.n.partCenter[1];
pos_z = curNode->u.n.partCenter[2];
}
/////////////////////////////////////////////////////////////////////////
//
// Process the active list window adding children to end of list
// Valid particles and accepted nodes are copied to the interact list
//
int hasChildren = 0;
int pcount = 0;
int ncount = 0;
for (int indx = begIndx; indx < endIndx; indx++) {
// If the current active element is a cell it will be
// ACCEPTED and copied to the interact list
// OPENED and its children will be added to the end of the active list
// IGNORED because it is too far away
if ((*curActive)[indx] >= this->nodeOffset) {
hasChildren = 1;
FNode* actNode = &this->fNode[(*curActive)[indx] - this->nodeOffset];
partRadius = actNode->u.n.partRadius;
distToNearPoint = distanceToNearestPoint(pos_x, pos_y, pos_z, actNode);
dx = actNode->u.n.partCenter[0] - pos_x;
dy = actNode->u.n.partCenter[1] - pos_y;
dz = actNode->u.n.partCenter[2] - pos_z;
r = sqrt(dx * dx + dy * dy + dz * dz);
// Node is ignored if it is too far away from the particle
// Distance from particle to particle radius exceeds critical radius
// Distance from particle to nearest side of node exceeds critical radius
if ((r - partRadius) > critRadius || distToNearPoint > critRadius) {
// Ignore node, move on to sibling of this node
}
else {
if (2*partRadius > (r * bhAngle)) {
// Open node, move on to first child
ID_T child =
this->fNode[(*curActive)[indx] - this->nodeOffset].u.n.nextNode;
while (child != -1) {
if (child >= this->nodeOffset) {
// Child is a node which is active and must be considered
newActive->push_back(child);
child = this->fNode[child - this->nodeOffset].u.n.sibling;
}
else {
// Child is a particle, add to interaction list
partInteract->push_back(child);
pcount++;
child = this->fParticle[child].sibling;
}
}
} else {
// Accept node, add to interact list, move on to sibling
nodeInteract->push_back((*curActive)[indx]);
ncount++;
}
}
}
}
/////////////////////////////////////////////////////////////////////////
//
// At this point a new level of children may have been added to active list
// Process children by dividing the node size and recursing
//
if (hasChildren) {
// Current item on active list is a cell
if (curId >= this->nodeOffset) {
// Process each child
ID_T child = fNode[curId - this->nodeOffset].u.n.nextNode;
while (child != -1) {
// Child is a node
if (child >= this->nodeOffset) {
FNode* childNode = &this->fNode[child - this->nodeOffset];
// Recurse on walk tree to process child
walkTreeBarnesAdjust(newActive, partInteract, nodeInteract,
child, bhAngle, critRadius);
child = childNode->u.n.sibling;
}
// Child is a particle
else {
walkTreeBarnesAdjust(newActive, partInteract, nodeInteract,
child, bhAngle, critRadius);
child = this->fParticle[child].sibling;
}
}
}
// Current item on active list is a particle
else {
walkTreeBarnesAdjust(newActive, partInteract, nodeInteract,
curId, bhAngle, critRadius);
}
}
/////////////////////////////////////////////////////////////////////////
//
// If no new items were added to active list we are done and can process
// the interact list for this particle p (which can't be a cell)
//
else {
if (curId > this->nodeOffset)
cout << "ERROR: POP OUT ON NODE " << curId << endl;
// Since the interact lists might contain accepted nodes from upper levels
// which need to be opened for this particle, adjust the lists first
vector<ID_T>* adjNodeInteract = new vector<ID_T>;
vector<ID_T>* adjPartInteract = new vector<ID_T>;
static Timings::TimerRef adjtimer = Timings::getTimer("Barnes Adjustment");
Timings::startTimer(adjtimer);
adjustInteraction(curId,
partInteract, nodeInteract,
adjPartInteract, adjNodeInteract,
bhAngle, critRadius);
Timings::stopTimer(adjtimer);
// Calculate force for the particle
this->fParticle[curId].force =
forceCalculation(curId, adjPartInteract, adjNodeInteract);
delete adjNodeInteract;
delete adjPartInteract;
}
// Active list is new for every recursion level
// Interact lists are appended to at each recursion level
// So interact lists must be popped by the correct number for this recursion
for (int i = 0; i < pcount; i++)
partInteract->pop_back();
for (int i = 0; i < ncount; i++)
nodeInteract->pop_back();
delete newActive;
}
/////////////////////////////////////////////////////////////////////////
//
// Recursion enters with a guess for the interact lists which were set
// on previous invocations of the method. Check the node interact list
// which contains nodes which were accepted to see if they should
// actually be opened relative to this new current particle or node
// If so remove from the nodeInteract list and add to the active list
//
// Particles in the interact list might actually be grouped and used
// with their parent node as an accept, but leaving them will lead to
// a better answer, not a worse. So we won't change the partInteract
//
/////////////////////////////////////////////////////////////////////////
void BHForceTree::adjustInteraction(
ID_T p0,
vector<ID_T>* partInteract,
vector<ID_T>* nodeInteract,
vector<ID_T>* adjPartInteract,
vector<ID_T>* adjNodeInteract,
POSVEL_T bhAngle,
POSVEL_T critRadius)
{
POSVEL_T dx, dy, dz, r, partRadius, distToNearPoint;
// Get location of particle being adjusted for
POSVEL_T pos_x = this->xx[p0];
POSVEL_T pos_y = this->yy[p0];
POSVEL_T pos_z = this->zz[p0];
// Copy all particles to the adjust list, will only add new particles
int numberOfParticles = (int) partInteract->size();
for (int p = 0; p < numberOfParticles; p++)
adjPartInteract->push_back((*partInteract)[p]);
// Process each node to see if status changes from accept to ignore or open
int numberOfNodes = (int) nodeInteract->size();
for (int n = 0; n < numberOfNodes; n++) {
FNode* curNode = &this->fNode[(*nodeInteract)[n] - this->nodeOffset];
partRadius = curNode->u.n.partRadius;
distToNearPoint = distanceToNearestPoint(pos_x, pos_y, pos_z, curNode);
dx = curNode->u.n.partCenter[0] - pos_x;
dy = curNode->u.n.partCenter[1] - pos_y;
dz = curNode->u.n.partCenter[2] - pos_z;
r = sqrt(dx * dx + dy * dy + dz * dz);
// Node is ignored if it is too far away from the particle
// Distance from particle to particle radius exceeds critical radius
// Distance from particle to nearest side of node exceeds critical radius
if ((r - partRadius) > critRadius || distToNearPoint > critRadius) {
// Ignore node, move on to sibling of this node
}
else {
if (2*partRadius > (r * bhAngle)) {
// Node must be opened and constituent parts examined
adjustNodeInteract(p0, curNode, adjPartInteract, adjNodeInteract,
bhAngle, critRadius);
} else {
// Accept node, add to interact list, move on to sibling
adjNodeInteract->push_back((*nodeInteract)[n]);
}
}
}
}
/////////////////////////////////////////////////////////////////////////
//
// Recursive part of interaction adjustment
// Examine children of current node recursively for inclusion into interaction
//
/////////////////////////////////////////////////////////////////////////
void BHForceTree::adjustNodeInteract(
ID_T p0,
FNode* curNode,
vector<ID_T>* adjPartInteract,
vector<ID_T>* adjNodeInteract,
POSVEL_T bhAngle,
POSVEL_T critRadius)
{
POSVEL_T dx, dy, dz, r, partRadius, distToNearPoint;
// Get location of particle being adjusted for
POSVEL_T pos_x = this->xx[p0];
POSVEL_T pos_y = this->yy[p0];
POSVEL_T pos_z = this->zz[p0];
// Current node is to be opened and recursively checked for interactions
ID_T child = curNode->u.n.nextNode;
while (child != -1) {
if (child < this->nodeOffset) {
// Child is a particle, add to adjusted particle interact list
adjPartInteract->push_back(child);
child = this->fParticle[child].sibling;
}
else {
// Child is a node, check to see if it should be opened, accepted, ignored
FNode* childNode = &this->fNode[child - this->nodeOffset];
partRadius = childNode->u.n.partRadius;
distToNearPoint = distanceToNearestPoint(pos_x, pos_y, pos_z, childNode);
dx = childNode->u.n.partCenter[0] - pos_x;
dy = childNode->u.n.partCenter[1] - pos_y;
dz = childNode->u.n.partCenter[2] - pos_z;
r = sqrt(dx * dx + dy * dy + dz * dz);
// Node is ignored if it is too far away from the particle
// Distance from particle to particle radius exceeds critical radius
// Distance from particle to nearest side of node exceeds critical radius
if ((r - partRadius) > critRadius || distToNearPoint > critRadius) {
// Ignore node, move on to sibling of this node
}
else {
if (2*partRadius > (r * bhAngle)) {
// Node must be opened and constituent parts examined
adjustNodeInteract(p0, childNode, adjPartInteract, adjNodeInteract,
bhAngle, critRadius);
} else {
// Accept node
adjNodeInteract->push_back(child);
}
}
child = this->fNode[child - this->nodeOffset].u.n.sibling;
}
}
}
/////////////////////////////////////////////////////////////////////////
//
// Short range gravity calculation for every particle in the tree
// Recurses through the tree saving previous work for reuse when popping
// out of recursion. Based on Barnes treecode with quick scan.
//
/////////////////////////////////////////////////////////////////////////
void BHForceTree::treeForceBarnesQuick(
POSVEL_T bhAngle, // Open a node
POSVEL_T critRadius) // Accept or ignore node not opened
{
ID_T root = this->particleCount;
vector<ID_T>* active = new vector<ID_T>;
vector<ID_T>* partInteract = new vector<ID_T>;
vector<ID_T>* nodeInteract = new vector<ID_T>;
active->push_back(root);
// Quick walk of tree accepts nodes that do not touch target node
walkTreeBarnesQuick(active, partInteract, nodeInteract,
root, bhAngle, critRadius);
delete active;
delete partInteract;
delete nodeInteract;
}
/////////////////////////////////////////////////////////////////////////
//
// Walk the BH tree for the given particle or node (identifier curId)
// Recursion starts with a new active list which will contain particles
// and nodes which possibly will contribute to the force on a particle.
// Particles on the active list will always be chosen for the interact list.
// Nodes on the active list may be OPENED if they are close enough
// or ACCEPTED and used in summary if they are not. Nodes that are opened
// have all their children (particles or nodes) added to the active list.
//
// After the children are added a new level of recursion starts by
// calculating a new size for that level, starting a fresh active list
// and building on the current interact lists.
//
// Recursion continues until the active list has been completely processed.
// When a level of recursion is complete the active list is destroyed
// and new items put on the interact lists are popped off.
//
// The advantage to this method is that items in the interaction list may
// not need to be processed again when we are doing the low levels of the tree.
//
/////////////////////////////////////////////////////////////////////////
void BHForceTree::walkTreeBarnesQuick(
vector<ID_T>* curActive, // nodes to be acted on
vector<ID_T>* partInteract, // particles for force
vector<ID_T>* nodeInteract, // nodes for force
ID_T curId, // current particle or node
POSVEL_T bhAngle, // open node
POSVEL_T critRadius) // accept or ignore
{
POSVEL_T dx, dy, dz, r, partRadius, distToNearPoint;
// Current active list
int begIndx = 0;
int endIndx = curActive->size();
// Construct active list for each recursion
vector<ID_T>* newActive = new vector<ID_T>;
// Set the location for the particle or node for the walk
POSVEL_T pos_x, pos_y, pos_z;
if (curId < this->nodeOffset) {
pos_x = this->xx[curId];
pos_y = this->yy[curId];
pos_z = this->zz[curId];
} else {
FNode* curNode = &this->fNode[curId - this->nodeOffset];
pos_x = curNode->u.n.partCenter[0];
pos_y = curNode->u.n.partCenter[1];
pos_z = curNode->u.n.partCenter[2];
}
/////////////////////////////////////////////////////////////////////////
//
// Process the active list window adding children to end of list
// Valid particles and accepted nodes are copied to the interact list
//
int hasChildren = 0;
int pcount = 0;
int ncount = 0;
for (int indx = begIndx; indx < endIndx; indx++) {
// If the current active element is a cell it will be
// ACCEPTED and copied to the interact list
// OPENED and its children will be added to the end of the active list
// IGNORED because it is too far away
if ((*curActive)[indx] >= this->nodeOffset) {
hasChildren = 1;
FNode* actNode = &this->fNode[(*curActive)[indx] - this->nodeOffset];
partRadius = actNode->u.n.partRadius;
distToNearPoint = distanceToNearestPoint(pos_x, pos_y, pos_z, actNode);
dx = actNode->u.n.partCenter[0] - pos_x;
dy = actNode->u.n.partCenter[1] - pos_y;
dz = actNode->u.n.partCenter[2] - pos_z;
r = sqrt(dx * dx + dy * dy + dz * dz);
// Node is ignored if it is too far away from the particle
// Distance from particle to particle radius exceeds critical radius
// Distance from particle to nearest side of node exceeds critical radius
if ((r - partRadius) > critRadius || distToNearPoint > critRadius) {
// Ignore node, move on to sibling of this node
}
else {
if (2*partRadius > (r * bhAngle)) {
// Open node, move on to first child
ID_T child =
this->fNode[(*curActive)[indx] - this->nodeOffset].u.n.nextNode;
while (child != -1) {
if (child >= this->nodeOffset) {
// Child is a node which is active and must be considered
newActive->push_back(child);
child = this->fNode[child - this->nodeOffset].u.n.sibling;
}
else {
// Child is a particle, add to interaction list
partInteract->push_back(child);
pcount++;
child = this->fParticle[child].sibling;
}
}
} else {
// Accept node, add to interact list, move on to sibling
nodeInteract->push_back((*curActive)[indx]);
ncount++;
}
}
}
}
/////////////////////////////////////////////////////////////////////////
//
// At this point a new level of children may have been added to active list
// Process children by dividing the node size and recursing
//
if (hasChildren) {
// Current item on active list is a cell
if (curId >= this->nodeOffset) {
// Process each child
ID_T child = fNode[curId - this->nodeOffset].u.n.nextNode;
while (child != -1) {
// Child is a node
if (child >= this->nodeOffset) {
FNode* childNode = &this->fNode[child - this->nodeOffset];
// Recurse on walk tree to process child
walkTreeBarnesQuick(newActive, partInteract, nodeInteract,
child, bhAngle, critRadius);
child = childNode->u.n.sibling;
}
// Child is a particle
else {
walkTreeBarnesQuick(newActive, partInteract, nodeInteract,
child, bhAngle, critRadius);
child = this->fParticle[child].sibling;
}
}
}
// Current item on active list is a particle
else {
walkTreeBarnesQuick(newActive, partInteract, nodeInteract,
curId, bhAngle, critRadius);
}
}
/////////////////////////////////////////////////////////////////////////
//
// If no new items were added to active list we are done and can process
// the interact list for this particle p (which can't be a cell)
//
else {
if (curId > this->nodeOffset)
cout << "ERROR: POP OUT ON NODE " << curId << endl;
this->fParticle[curId].force =
forceCalculation(curId, partInteract, nodeInteract);
}
// Active list is new for every recursion level
// Interact lists are appended to at each recursion level
// So interact lists must be popped by the correct number for this recursion
for (int i = 0; i < pcount; i++)
partInteract->pop_back();
for (int i = 0; i < ncount; i++)
nodeInteract->pop_back();
delete newActive;
}
/////////////////////////////////////////////////////////////////////////
//
// Short range force calculation
// Potential is calculated and is used to determine the acceleration of
// the particle. Acceleration is applied to the current velocity to
// produce the velocity at the next time step.
//
/////////////////////////////////////////////////////////////////////////
POSVEL_T BHForceTree::forceCalculation(
ID_T p0, // Target particle index
vector<ID_T>* partInteract, // Particles acting on p
vector<ID_T>* nodeInteract) // Nodes acting on p
{
POSVEL_T accel[DIMENSION];
POSVEL_T phi = 0.0;
POSVEL_T pos0_x = this->xx[p0];
POSVEL_T pos0_y = this->yy[p0];
POSVEL_T pos0_z = this->zz[p0];
for (int dim = 0; dim < DIMENSION; dim++)
accel[dim] = 0.0;
int numberOfNodes = (int) nodeInteract->size();
int numberOfParticles = (int) partInteract->size();
// Particles contributing to the force use location and mass of one particle
for (int p = 0; p < numberOfParticles; p++) {
ID_T particle = (*partInteract)[p];
if (p0 != particle) {
POSVEL_T dx = this->xx[particle] - pos0_x;
POSVEL_T dy = this->yy[particle] - pos0_y;
POSVEL_T dz = this->zz[particle] - pos0_z;
POSVEL_T r2 = dx * dx + dy * dy + dz * dz;
POSVEL_T f_over_r = this->mass[particle] * m_fl->f_over_r(r2);
//POSVEL_T f_over_r = this->mass[particle] / r2;
phi -= f_over_r;
//if (p0 == 171893) cout << "Top Particle used " << particle << " phi " << phi << endl;
accel[0] += dx * f_over_r * m_fcoeff;
accel[1] += dy * f_over_r * m_fcoeff;
accel[2] += dz * f_over_r * m_fcoeff;
this->vx[p0] += dx * f_over_r * m_fcoeff;
this->vy[p0] += dy * f_over_r * m_fcoeff;
this->vz[p0] += dz * f_over_r * m_fcoeff;
}
}
// Nodes contributing to force use center of mass and total particle mass
for (int n = 0; n < numberOfNodes; n++) {
FNode* node = &this->fNode[(*nodeInteract)[n] - this->nodeOffset];
POSVEL_T dx = node->u.n.partCenter[0] - pos0_x;
POSVEL_T dy = node->u.n.partCenter[1] - pos0_y;
POSVEL_T dz = node->u.n.partCenter[2] - pos0_z;
POSVEL_T r2 = dx * dx + dy * dy + dz * dz;
POSVEL_T f_over_r = node->u.n.partMass * m_fl->f_over_r(r2);
//POSVEL_T f_over_r = node->u.n.partMass / r2;
phi -= f_over_r;
//if (p0 == 171893) cout << "Top node used " << (*nodeInteract)[n] << " phi " << phi << endl;
accel[0] += dx * f_over_r * m_fcoeff;
accel[1] += dy * f_over_r * m_fcoeff;
accel[2] += dz * f_over_r * m_fcoeff;
this->vx[p0] += dx * f_over_r * m_fcoeff;
this->vy[p0] += dy * f_over_r * m_fcoeff;
this->vz[p0] += dz * f_over_r * m_fcoeff;
}
return phi;
}
POSVEL_T BHForceTree::forceCalculationFast(
ID_T p0, // Target particle index
vector<POSVEL_T>* xInteract,
vector<POSVEL_T>* yInteract,
vector<POSVEL_T>* zInteract,
vector<POSVEL_T>* mInteract)
{
POSVEL_T phi = 0.0;
POSVEL_T pos0_x = this->xx[p0];
POSVEL_T pos0_y = this->yy[p0];
POSVEL_T pos0_z = this->zz[p0];
int nInteract = (int) xInteract->size();
for (int p = 0; p < nInteract; p++) {
POSVEL_T dx = (*xInteract)[p] - pos0_x;
POSVEL_T dy = (*yInteract)[p] - pos0_y;
POSVEL_T dz = (*zInteract)[p] - pos0_z;
POSVEL_T r2 = dx * dx + dy * dy + dz * dz;
POSVEL_T f_over_r = (*mInteract)[p] * m_fl->f_over_r(r2);
//POSVEL_T f_over_r = this->mass[particle] / r2;
this->vx[p0] += dx * f_over_r * m_fcoeff;
this->vy[p0] += dy * f_over_r * m_fcoeff;
this->vz[p0] += dz * f_over_r * m_fcoeff;
}
return phi;
}
/////////////////////////////////////////////////////////////////////////
//
// Return the distance^2 from location to the closest point on FNode
//
/////////////////////////////////////////////////////////////////////////
POSVEL_T BHForceTree::distanceToNearestPoint(
POSVEL_T pos_x,
POSVEL_T pos_y,
POSVEL_T pos_z,
FNode* node)
{
// Calculate bounding box of current node
// Nearest point in bounding box decides whether particle or node is used
POSVEL_T dx, dy, dz, r;
POSVEL_T minBound[DIMENSION], maxBound[DIMENSION];
for (int dim = 0; dim < DIMENSION; dim++) {
minBound[dim] = node->geoCenter[dim] - (node->geoSide[dim] * 0.5);
maxBound[dim] = node->geoCenter[dim] + (node->geoSide[dim] * 0.5);
}
if (pos_x < minBound[0])
dx = minBound[0] - pos_x;
else if (pos_x > maxBound[0])
dx = pos_x - maxBound[0];
else
dx = 0.0;
if (pos_y < minBound[1])
dy = minBound[1] - pos_y;
else if (pos_y > maxBound[1])
dy = pos_y - maxBound[1];
else
dy = 0.0;
if (pos_z < minBound[2])
dz = minBound[2] - pos_z;
else if (pos_z > maxBound[2])
dz = pos_z - maxBound[2];
else
dz = 0.0;
r = sqrt(dx * dx + dy * dy + dz * dz);
return r;
}
/////////////////////////////////////////////////////////////////////////
//
// Return the distance from location to the fNode center of mass
//
/////////////////////////////////////////////////////////////////////////
POSVEL_T BHForceTree::distanceToCenterOfMass(
POSVEL_T xLoc,
POSVEL_T yLoc,
POSVEL_T zLoc,
FNode* node)
{
POSVEL_T xdist = (POSVEL_T) fabs(xLoc - node->u.n.partCenter[0]);
POSVEL_T ydist = (POSVEL_T) fabs(yLoc - node->u.n.partCenter[1]);
POSVEL_T zdist = (POSVEL_T) fabs(zLoc - node->u.n.partCenter[2]);
POSVEL_T dist = sqrt((xdist * xdist) + (ydist * ydist) + (zdist * zdist));
return dist;
}
/////////////////////////////////////////////////////////////////////////
//
// Return the distance from location to the fNode furthest corner
//
/////////////////////////////////////////////////////////////////////////
POSVEL_T BHForceTree::distanceToFarCorner(
POSVEL_T xLoc,
POSVEL_T yLoc,
POSVEL_T zLoc,
FNode* node)
{
POSVEL_T distance = 0.0;
POSVEL_T corner[DIMENSION];
POSVEL_T xdist, ydist, zdist, dist;
for (int k = -1; k <= 1; k=k+2) {
corner[2] = node->geoCenter[2] + (k * (node->geoSide[2] * 0.5));
for (int j = -1; j <= 1; j=j+2) {
corner[1] = node->geoCenter[1] + (j * (node->geoSide[1] * 0.5));
for (int i = -1; i <= 1; i=i+2) {
corner[0] = node->geoCenter[0] + (i * (node->geoSide[0] * 0.5));
xdist = (POSVEL_T) fabs(xLoc - corner[0]);
ydist = (POSVEL_T) fabs(yLoc - corner[1]);
zdist = (POSVEL_T) fabs(zLoc - corner[2]);
dist = sqrt((xdist * xdist) + (ydist * ydist) + (zdist * zdist));
if (dist > distance)
distance = dist;
}
}
}
return distance;
}
/////////////////////////////////////////////////////////////////////////
//
// Return the distance from location to the fNode nearest corner
//
/////////////////////////////////////////////////////////////////////////
POSVEL_T BHForceTree::distanceToNearCorner(
POSVEL_T xLoc,
POSVEL_T yLoc,
POSVEL_T zLoc,
FNode* node)
{
POSVEL_T distance = MAX_FLOAT;
POSVEL_T corner[DIMENSION];
POSVEL_T xdist, ydist, zdist, dist;
for (int k = -1; k <= 1; k=k+2) {
corner[2] = node->geoCenter[2] + (k * (node->geoSide[2] * 0.5));
for (int j = -1; j <= 1; j=j+2) {
corner[1] = node->geoCenter[1] + (j * (node->geoSide[1] * 0.5));
for (int i = -1; i <= 1; i=i+2) {
corner[0] = node->geoCenter[0] + (i * (node->geoSide[0] * 0.5));
xdist = (POSVEL_T) fabs(xLoc - corner[0]);
ydist = (POSVEL_T) fabs(yLoc - corner[1]);
zdist = (POSVEL_T) fabs(zLoc - corner[2]);
dist = sqrt((xdist * xdist) + (ydist * ydist) + (zdist * zdist));
if (dist < distance)
distance = dist;
}
}
}
return distance;
}
/////////////////////////////////////////////////////////////////////////
//
// Print BH tree with indentations indicating levels
// Since the tree has been threaded changing the recursive tree with children
// into an iterative tree with next nodes and parents, walk the tree
// iteratively keeping track of parents to indicate when levels change
//
/////////////////////////////////////////////////////////////////////////
void BHForceTree::printBHForceTree()
{
ID_T curIndex = this->nodeOffset;
vector<ID_T> parents;
parents.push_back(-1);
ID_T parentIndex = 0;
while (curIndex != -1) {
// Get the parent of the current index
ID_T parent;
if (curIndex >= this->nodeOffset)
parent = this->fNode[curIndex - this->nodeOffset].u.n.parent;
else
parent = this->fParticle[curIndex].parent;
// Pop the stack of parents until the level is right
while (parent != parents[parentIndex]) {
parents.pop_back();
parentIndex--;
}
// Print FNode
if (curIndex >= this->nodeOffset) {
FNode* fn = &this->fNode[curIndex-this->nodeOffset];
cout << parentIndex << ":" << setw(parentIndex) << " ";
cout << "N " << curIndex
<< " sibling " << fn->u.n.sibling
<< " next " << fn->u.n.nextNode
<< " parent " << fn->u.n.parent
<< " [" << (fn->geoCenter[0]-fn->geoSide[0]/2.0)
<< ":" << (fn->geoCenter[0]+fn->geoSide[0]/2.0) << "] "
<< " [" << (fn->geoCenter[1]-fn->geoSide[1]/2.0)
<< ":" << (fn->geoCenter[1]+fn->geoSide[1]/2.0) << "] "
<< " [" << (fn->geoCenter[2]-fn->geoSide[2]/2.0)
<< ":" << (fn->geoCenter[2]+fn->geoSide[2]/2.0) << "] "
<< " (" << fn->u.n.partCenter[0]
<< " ," << fn->u.n.partCenter[1]
<< " ," << fn->u.n.partCenter[2]
<< ") MASS " << fn->u.n.partMass
<< " RADIUS " << fn->u.n.partRadius
<< endl;
// Push back the new FNode which will have children
parents.push_back(curIndex);
parentIndex++;
// Walk to next node (either particle or node)
curIndex = this->fNode[curIndex-this->nodeOffset].u.n.nextNode;
}
// Print FParticle
else {
cout << parentIndex << ":" << setw(parentIndex) << " ";
cout << "P " << curIndex
<< " sibling " << this->fParticle[curIndex].sibling
<< " next " << this->fParticle[curIndex].nextNode
<< " parent " << this->fParticle[curIndex].parent
<< " (" << xx[curIndex]
<< " ," << yy[curIndex]
<< " ," << zz[curIndex] << ")" << endl;
// Walk to next node (either particle or node)
curIndex = this->fParticle[curIndex].nextNode;
}
}
}
/////////////////////////////////////////////////////////////////////////
//
// Print the force values for comparison
//
/////////////////////////////////////////////////////////////////////////
void BHForceTree::printForceValues()
{
for (int p = 0; p < this->particleCount; p++) {
cout << "Particle: " << setw(8) << p
<< " force " << this->fParticle[p].force << endl;
}
}
/////////////////////////////////////////////////////////////////////////
//
// Get the index of the child which should contain this particle
//
/////////////////////////////////////////////////////////////////////////
int BHForceTree::getChildIndex(FNode* node, ID_T pindx)
{
// Vary Z dimension fastest in making octtree children
int index = 0;
if (this->xx[pindx] >= node->geoCenter[0]) index += 4;
if (this->yy[pindx] >= node->geoCenter[1]) index += 2;
if (this->zz[pindx] >= node->geoCenter[2]) index += 1;
return index;
}