我正在模拟两个星团的碰撞，所以基本上在C ++中做了一个N体问题。我正在使用Barnes-Hut方法来帮助实现这一目标。我编写了我的QuadTree代码和运行模拟的主要代码，直到201-250步骤，程序变得没有响应。我不确定这只是一个性能问题，还是代码中发生了令人讨厌的事情。有人可以查看它并解释它冻结的原因，或者如何提高性能？
以下是使用Barnes-Hut进行QuadTree构造和计算的代码：
#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);
}
};
请注意，这使用了一个名为Vector3D的类，它基本上用于物理矢量定义和操作。这是时间演变主循环的代码：
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;
}
这使用了执行Runge-Kutta计算的外部函数。这已被用于其他程序，所以不要担心VecFRK4功能。
可以找到运行此代码所需的完整代码集合here。
我能够修复并完成此代码。我将为任何可能在创建QuadTree时遇到类似问题的人发布我的解决方案。
问题在于QuadTree代码本身。在评论中，有人提到is_internal（）方法将始终返回false，因为在if语句创建新的子节点之前创建了具有0质量的新子节点，因此永远不会满足条件为真。下一个问题是无条件地创建新的子节点。我们只想在insert（）方法中第一次调用它们时创建子节点，我们只想在子节点不存在的情况下创建子节点（即我们不想覆盖已存在的子节点）。为了解决这个问题，我们认识到第一次使用子节点是在is_internal（）方法中。这取决于现有的子节点，因此我们必须在该方法中创建子节点：
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
}
现在，为了防止每次调用createSubnodes（）时覆盖新的子节点，我们只需在创建它们之前询问它们是否存在：
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
}
这样可以防止程序冻结，因为四叉树不会在没有约束的情况下完成。
下一个问题是这会导致内存泄漏。该程序将继续创建子QuadTrees并将其保留在内存中。我们需要创建一个删除方法，通过从底部开始删除子节点来释放空间。为此，我们可以：
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
}
在主程序中，在每个时间步的末尾调用此方法。另外，在每个时间步之后调用forces.clear()以删除不再需要的力并为新力量腾出空间，从而防止另一个内存泄漏。
我希望有一天有人觉得这很有帮助！