-1

I'm simulating the collision of two star clusters, so basically doing an N-body problem, in C++. I'm using the Barnes-Hut method to help do this. I've written my QuadTree code and the main code to run the simulation and it works up until the 201-250th step, where the program then becomes unresponsive. I am not sure if this is just a performance issue or if there's something nasty happening with the code. Can someone look over it and explain why it freezes, or how to improve the performance?

Here's the code for the QuadTree construction and calculation using Barnes-Hut:

#define threshold 1.0   // threshold for determining "far-ness"
                        // 0.5 for best results, higher means faster 
#define G 1.57633e-17   // ly^3 Msun^-1 s^-2, gravitational constant

struct Body {
    double mass;
    Vector3D position;
    Vector3D velocity;

    // makes a new body
    Body(double m, Vector3D pos, Vector3D vel) : mass(m), position(pos), velocity(vel) {};

    Body() : mass(0), position(0, 0, 0), velocity(0, 0, 0) {};
};

struct QuadTreeNode;

struct QuadTreeNode {

/* ---------- Stored Properties of this node ----------- */

    // the total mass and CoM for all bodies inserted
    // below this node in the tree
    Body body;

    // rescale positions to fit in nodes
    // note: will need to have all bodies in
    // at least positive coordinates to be included
    // body.position.update(body.position.Getx()/scale, 
    //                      body.position.Gety()/scale);
    // scale defined in program

    // the four subnodes / quadrants
    QuadTreeNode* A = NULL;
    QuadTreeNode* B = NULL;
    QuadTreeNode* C = NULL;
    QuadTreeNode* D = NULL;

    // location+size parameters of each node (center)
    double length;          // = 1.0 for depth = 0
    double horiz_offset;    // = 0.5*length for depth = 0
    double vert_offset;     // "                        "

    // used to initialize root node parameters
    bool root_node = true;

/* ----------------------------------------------------- */

    // search which subnode a body is in
    int Subnode(Body temp){
        double x = temp.position.Getx();
        double y = temp.position.Gety();

        // define subnode boundaries
        double node_xmax = this->horiz_offset + 0.5*this->length;
        double node_xmin = this->horiz_offset - 0.5*this->length;
        double node_ymax = this->vert_offset + 0.5*this->length;
        double node_ymin = this->vert_offset - 0.5*this->length;

        // 1 = first quadrant (A), 2 = second quadrant (B), etc.
        if       ((x > (node_xmax/2.)) && (x < node_xmax) && 
                  (y > (node_ymax/2.)) && (y < node_ymax)){return(1);
        }else if ((x > node_xmin) && (x < (node_xmax/2.)) && 
              (y > (node_ymax/2.)) && (y < node_ymax)){return(2);
        }else if ((x > node_xmin) && (x < (node_xmax/2.)) && 
              (y > node_ymin) && (y < (node_ymax/2.))){return(3);
        }else if ((x > (node_xmax/2.)) && (x < node_xmax) && 
              (y > node_ymin) && (y < (node_ymax/2.))){return(4);
        }else{
            temp.mass = 0.0; 
            return(0);
        }
    }

    // define new parameters of subnodes
    void createSubnodes(){
        // 1st quadrant / subnode
        this->A = new QuadTreeNode();
        this->A->length = this->length/2.;
        this->A->horiz_offset = this->horiz_offset + 0.25*this->length;
        this->A->vert_offset = this->vert_offset + 0.25*this->length;

        // 2nd quadrant / subnode
        this->B = new QuadTreeNode();
        this->B->length = this->length/2.;
        this->B->horiz_offset = this->horiz_offset - 0.25*this->length;
        this->B->vert_offset = this->vert_offset + 0.25*this->length;

        // 3rd quadrant / subnode
        this->C = new QuadTreeNode();
        this->C->length = this->length/2.;
        this->C->horiz_offset = this->horiz_offset - 0.25*this->length;
        this->C->vert_offset = this->vert_offset - 0.25*this->length;

        // 4th quadrant / subnode
        this->D = new QuadTreeNode();
        this->D->length = this->length/2.;
        this->D->horiz_offset = this->horiz_offset + 0.25*this->length;
        this->D->vert_offset = this->vert_offset - 0.25*this->length;
    }

    // insert body into specific subnode
    void insert_Subnode(Body temp){

        createSubnodes();

        // subnodes are not a root node
        this->A->root_node = false;
        this->B->root_node = false;
        this->C->root_node = false;
        this->D->root_node = false;

        // recursively insert body into new node 
        if (Subnode(temp) == 1) {
            this->A->insert(temp);
        }
        else if (Subnode(temp) == 2) {
            this->B->insert(temp);
        }
        else if (Subnode(temp) == 3) {
            this->C->insert(temp);
        }
        else if (Subnode(temp) == 4){
            this->D->insert(temp);
        }else return;
    }

    // check if current node is internal or external
    bool is_internal(){
        this->A = new QuadTreeNode();
        this->B = new QuadTreeNode();
        this->C = new QuadTreeNode();
        this->D = new QuadTreeNode();
        if (this->A->body.mass > 0) return true;
        if (this->B->body.mass > 0) return true;
        if (this->C->body.mass > 0) return true;
        if (this->D->body.mass > 0) return true;
        return false;   // this is an external node
    }

    // this is the main function for constructing the QuadTree
    // insert body into this node
    void insert(Body next) {    // consider the next body 
                                // to insert into this node

        // if there's no body in this node, put the next body in here
        if (this->body.mass == 0) {
            this->body = next;
            return;
        }

        // initialize to max values if this is a root node
        if(this->root_node){
            this->length = 1.0;
            this->horiz_offset = 0.5*this->length;
            this->vert_offset = 0.5*this->length;
        }

        // if this is an internal node
        if (is_internal()) {
            // update the node's body.location (CoM)
            this->body.position = (this->body.mass*this->body.position 
                       + next.mass*next.position)*(1./(this->body.mass + next.mass));
            // update the node's body.mass
            this->body.mass += next.mass;

            // recursively insert the next body into appropriate subnode
            insert_Subnode(next);           
        } else {    // this is an external node
                    // since this->body.mass is not 0, then
                    // this node contains a body
                    // so then recursively insert this body
                    // and the next body into their appropriate subnodes
            insert_Subnode(this->body);
            insert_Subnode(next);

            // then update this node's CoM and total mass
            this->body.position = (this->body.mass*this->body.position 
                      + next.mass*next.position)*(1./(this->body.mass + next.mass));
            // update the node's body.mass
            this->body.mass += next.mass;
        }

        return;

    }


// Now we calculate the force on a body

    // calculates gravitational force
    Vector3D Fg(Body temp){
        double M1, M2, dist;
        Vector3D  Dr, F;

        M1 = this->body.mass;
        M2 = temp.mass;

        // this takes the position vector difference
        Dr =  temp.position - this->body.position;
        // and this is the magnitude
        dist = Dr.GetMagnitude();

        if(dist == 0.0) return(Vector3D(0.0, 0.0, 0.0));    // do not divide by 0

        // here is the gravitational force
        F = (-G*M1*M2/(pow(dist,3.))) * Dr;

        return(F);
    }

    // net force on the next body
    Vector3D netforce_on(Body next) {

        Vector3D Fnet(0.0, 0.0, 0.0);
        Vector3D dvdt;

        if ((!is_internal()) && (this->body.mass != next.mass)){
            Fnet = Fnet + Fg(next);
        }
        else {
            Vector3D Dr;
            double dist;
            Dr =  next.position - this->body.position;
            dist = Dr.GetMagnitude();

            // parameter to decide what is 'far'
            double dist_ratio = (this->length) / dist;

            if (dist_ratio < threshold){ // far enough
                Fnet = Fnet + Fg(next);
            } else{ // too close, check subnode forces
                Fnet = Fnet + this->A->netforce_on(next);
                Fnet = Fnet + this->B->netforce_on(next);
                Fnet = Fnet + this->C->netforce_on(next);
                Fnet = Fnet + this->D->netforce_on(next);
            }
        }

        dvdt = (1./next.mass)*Fnet;

        return(dvdt);

    }

};

Note that this uses a class called Vector3D that is basically for physics vector definitions and operations. And here's the code for the main loop of the time evolution:

    int ii = 0;                     // index for printing every 50s
    for (double t = 0; t < Tmax; t += dt) {
        static int div = (int)(50/dt);      // print only every 50s

        // create a new QuadTree
        if(t==0) cerr << "Creating 1st QuadTree..." << endl;    // checking if this works the first time
//      QuadTreeNode tree;
        QuadTreeNode* tree;
        tree = new QuadTreeNode();  

        // insert all of the masses into the tree
        if(t==0) cerr << "Inserting Bodies into 1st QuadTree..." << endl;
        for (int i = 0; i < bodies.size(); i++) {
            if(t==0){
                if(i%50 == 0){
                    cerr << "Inserting " << i << "th body..." << endl;
                }
            }
            tree->insert(bodies[i]);
        }

        if(t==0) cerr << "1st Body insertion finished." << endl;
        if(t==0) cerr << "Beginning RK4 Calculation..." << endl;

        for(int j = 0; j < bodies.size(); j++){
            Vector3D rtold, vtold, rt, vt;

            if(ii%div == 0){
                outfile << t+dt << " ";
            }

            // compute force on each body
//          vector<Vector3D> forces;    // declared globally for f_v
            for (int i = 0; i < bodies.size(); i++) {
                forces.push_back(tree->netforce_on(bodies[i]));
            }

            /* These lines update the position and velocity of jth body */
            rtold = bodies[j].position;
            vtold = bodies[j].velocity;
            rt = rtold + VecFRK4(0,f_r,f_v,j,t,rtold,vtold,dt);
            vt = vtold + VecFRK4(1,f_r,f_v,j,t,rtold,vtold,dt);

            /* set the old values to the new updated values for next iteration */
            bodies[j].position = rt;
            bodies[j].velocity = vt;

            if(ii%div == 0){
                outfile << rt.Getx() << " " << rt.Gety() << " ";
            }

        }
        if(ii%div == 0){
        outfile << endl;
    }

        if(ii%(50) == 0){       // note every 1e4 seconds
            cerr << "Calculating... t = " << t+dt << " ..." << endl;
        }

        ii++;

        delete tree;

    }

This uses an external function that does a Runge-Kutta calculation. That's been used in other programs so don't worry about that VecFRK4 function.

The complete collection of code necessary to run this can be found here.

Zachary
  • 47
  • 5

1 Answers1

0

I was able to fix and complete this code. I will post my solution for any of those who may have similar problems creating a QuadTree.

The problem lies within the QuadTree code itself. In the comments, there was mention that the is_internal() method will always returns false like that since creating the new subnodes right before the if statements is creating fresh subnodes with 0 mass and thus the conditions to be true will never be satisfied. The next problem is creating new subnodes unconditionally. We only want to create the subnodes the first time they are called in the insert() method and we only want to create them if the subnodes are not already there (ie we don't want to overwrite already existing subnodes). To fix this, we recognize the first use of subnodes to be in the is_internal() method. This depends on existing subnodes so we must create the subnodes in that method:

bool is_internal(){

    createSubnodes();

    if (this->A->body.mass > 0) return true;
    if (this->B->body.mass > 0) return true;
    if (this->C->body.mass > 0) return true;
    if (this->D->body.mass > 0) return true;
    return false;   // this is an external node
}

Now, to prevent new subnodes being overwritten each time createSubnodes() is called, we simply ask if they exist before creating them:

void createSubnodes(){
    // 1st quadrant / subnode
    if (!this->A) {    // create and initialize subnode if there's
                       // not an existing one already
        this->A = new QuadTreeNode();
        this->A->length = this->length/2.;
        this->A->horiz_offset = this->horiz_offset + 0.25*this->length;
        this->A->vert_offset = this->vert_offset + 0.25*this->length;
    }
    // repeat for B, C, and D subnodes
}

This will prevent the freezing of the program since quadtrees won't just be made without bound now.

The next problem is that this will cause a memory leak. The program will keep creating sub QuadTrees and keeping them in the memory. We need to create a delete method to free up the space by deleting the subnodes starting from the bottom. To do this, we can make:

void free() {       // delete the quadtree after use

    if (this->A) {                  // if this subnode exists,
        this->A->free();            // delete it's subnodes (keeps checking)
        delete this->A;             // until a NULL subnode is reached
    }
    // repeat for B, C, D
}

and in the main program, call this method at the end of each time step. Additionally, call forces.clear() after each time step as well to delete the no longer needed forces and make room for the new forces, thus preventing another memory leak.

I hope someone finds this helpful one day!

Zachary
  • 47
  • 5