C++中的N体仿真

Question

我正在尝试实现2维n体仿真的OpenMP版本.

但是存在一个问题:我假设每个粒子的初始速度和加速度都为零.当颗粒首先聚集在一起时,它们会高速分散,不再聚集.

这似乎与牛顿法不一致,对吧？

有人可以解释为什么会发生这种情况以及如何解决错误？

这是我的代码的一部分:

/* update one frame */
void update() {
    int i, j;

    omp_set_num_threads(8);
    # pragma omp parallel private(j)
    {

    # pragma omp for schedule(static)
        for ( i = 0; i < g_N; ++i ) {
            g_pv[i].a_x = 0.0;
            g_pv[i].a_y = 0.0;
            for ( j = 0; j < g_N; ++j ) {
                if (i == j)
                    continue;
                double r_2 = pow((g_pv[i].pos_x - g_pv[j].pos_x),2) + pow((g_pv[i].pos_y - g_pv[j].pos_y),2);
                g_pv[i].a_x += (-1) * G * g_pv[j].m * (g_pv[i].pos_x - g_pv[j].pos_x) / (pow(r_2 + e,1.5));
                g_pv[i].a_y += (-1) * G * g_pv[j].m * (g_pv[i].pos_y - g_pv[j].pos_y) / (pow(r_2 + e,1.5));             
            } 
            g_pv[i].v_x += period * g_pv[i].a_x;
            g_pv[i].v_y += period * g_pv[i].a_y;
        }
    # pragma omp for schedule(static)
        for ( int i = 0; i < g_N; ++i ) {
            g_pv[i].pos_x += g_pv[i].v_x * period;
            g_pv[i].pos_y += g_pv[i].v_y * period;
        }
    }
}

不要担心OpenMP,只需将其视为顺序版本.OpenMP不会对结果产生太大影响.

编辑:澄清,这是整个代码(这部分可能有一些错误,但我描述的问题应该出现在上面的代码部分)

# include <iostream>
# include <fstream>
# include <iomanip>
# include <cmath>
# include <vector>
# include <cstdlib>
# include <omp.h>

# include <GL/glew.h>
# include <GL/freeglut.h>
# include <GL/gl.h>

using namespace std;

/* the size of the opengl window */
# define WIDTH 2000
# define HEIGHT 2000

/* define the global constants */
const double G = 6.67 * pow(10, -11);
// const double G = 6.67;
const double e = 0.00001;
const double period = 1;

/* define the structure of particle */
struct particle
{
    double m;
    double pos_x;
    double pos_y;
    double v_x;
    double v_y;
    double a_x;
    double a_y;

    particle(double m = 0, double pos_x = 0, double pos_y = 0, 
            double v_x = 0, double v_y = 0, double a_x = 0, double a_y = 0)
    {
        this->m         = m;
        this->pos_x     = pos_x;
        this->pos_y     = pos_y;
        this->v_x       = v_x;
        this->v_y       = v_y;
        this->a_x       = a_x;
        this->a_y       = a_y;
    }
};

/* define the global data */
int g_N;                        // number of particles
vector<particle> g_pv;          // particle vector

void setUp();
void update();
void display();

int main(int argc, char ** argv) {

    /* set up the window */
    glutInit(&argc, argv);
    glutInitDisplayMode (GLUT_SINGLE | GLUT_RGB);
    glutInitWindowSize (WIDTH, HEIGHT);
    glutInitWindowPosition (100, 100);
    glutCreateWindow ("openmp");

    /* initialize */
    setUp();

    glutDisplayFunc(display);
    glutMainLoop();

    return 0;
}

/* read the input data */
void setUp() {
    glMatrixMode (GL_PROJECTION);
    glLoadIdentity ();
    /* Sets a 2d projection matrix
     * (0,0) is the lower left corner (WIDTH, HEIGHT) is the upper right */
    glOrtho (0, WIDTH, 0, HEIGHT, 0, 1);
    glDisable(GL_DEPTH_TEST);

    ifstream inFile;
    inFile.open("input_25.txt");

    inFile >> g_N;
    g_pv.resize(g_N);
    for ( int i = 0; i < g_N; ++i )
    {
        inFile >> g_pv[i].m >> g_pv[i].pos_x >> g_pv[i].pos_y
               >> g_pv[i].v_x >> g_pv[i].v_y >> g_pv[i].a_x >> g_pv[i].a_y; 
    }
    inFile.close();
}

/* display in openGL */
void display(void) {
    glClear(GL_COLOR_BUFFER_BIT);
    for(int i = 0; i < g_pv.size(); ++i) {
        /* Get the ith particle */
        particle p = g_pv[i];

        /* Draw the particle as a little square. */
        glBegin(GL_QUADS);
        glColor3f (1.0, 1.0, 1.0);
        glVertex2f(p.pos_x + 2, p.pos_y + 2);
        glVertex2f(p.pos_x - 2, p.pos_y + 2);
        glVertex2f(p.pos_x - 2, p.pos_y - 2);
        glVertex2f(p.pos_x + 2, p.pos_y - 2);
        glEnd();
    }   
    update();
    glutPostRedisplay();
    glFlush ();
}

/* update one frame */
void update() {
    int i, j;

    omp_set_num_threads(8);
    # pragma omp parallel private(j)
    {
        /* compute the force */
        # pragma omp for schedule(static)
            for ( i = 0; i < g_N; ++i ) {
                g_pv[i].a_x = 0.0;
                g_pv[i].a_y = 0.0;
                for ( j = 0; j < g_N; ++j ) {
                    if (i == j)
                        continue;
                    double r_2 = pow((g_pv[i].pos_x - g_pv[j].pos_x),2) + pow((g_pv[i].pos_y - g_pv[j].pos_y),2);
                    g_pv[i].a_x += (-1) * G * g_pv[j].m * (g_pv[i].pos_x - g_pv[j].pos_x) / (pow(r_2 + e,1.5));
                    g_pv[i].a_y += (-1) * G * g_pv[j].m * (g_pv[i].pos_y - g_pv[j].pos_y) / (pow(r_2 + e,1.5));             
                } 
                g_pv[i].v_x += period * g_pv[i].a_x;
                g_pv[i].v_y += period * g_pv[i].a_y;
            }

        /* compute the velocity */
        # pragma omp for schedule(static)
            for ( int i = 0; i < g_N; ++i ) {
                g_pv[i].pos_x += g_pv[i].v_x * period;
                g_pv[i].pos_y += g_pv[i].v_y * period;
            }
    }
}

Answer 1

我正在将我的评论扩展到一个答案(正如Z boson所建议的那样),并提出了一些关于如何解决问题的建议.

这个问题真的属于Computational Science.SE,因为我认为代码本身没有任何问题,但算法似乎有问题:随着粒子越来越近,你可以得到一个力量G / pow(e,1.5) ~ G * 1e7.这很大.非常非常大(与您的时间步长相比).为什么？假设你有两个行星,一个坐在(0,0),一个小坐在(0,1),后者对前者的加速度非常大.在下一步中,小行星将处于(0,-100)或其他任何位置,并且其上的力为零,这意味着它将永远不会返回,现在也具有相当大的速度.你的模拟已经爆炸了.

如你所说,这与牛顿定律并不相符,所以这表明你的数值方案已经失败了.不用担心,有几种方法可以解决这个问题.您已经预料到了,就像您添加的那样e.比如说10,让它变大,你应该没事.或者,设置一个非常小的时间步长period.如果粒子太接近或者如果行星碰撞会发生什么(可能是爆炸和消失,或者抛出异常),你也可以不计算力量.或者说具有r ^{- 2}势或g_pv[i].a_y += (-1) * G * g_pv[j].m * (g_pv[i].pos_y - g_pv[j].pos_y) * (1 / pow(r_2 + e,1.5) - 1 / pow(r_2 + e,2.5)); 类似物的排斥力: 这类似于现象学的Lennard-Jones相互作用结合了由泡利不相容原理引起的排斥.请注意,您可以增加排斥的锐度(如果您愿意,可以2.5改为12.5),这意味着排斥效果较差(好),但需要更准确地解决,导致更小period(不好) ).理想情况下,您将从不会导致冲突的初始配置开始,但这几乎是不可能预测的.最后,您可能希望使用上面列出的方法的组合.

使用OpenMP可能会导致轻微的加速,但您应该使用算法来实现远程力,例如Barnes-Hut模拟.有关更多信息,请参阅2013年暑期学校关于复杂粒子系统中远程相互作用的快速方法以及他们在最新发展中发布的小册子(免费提供).您也可能不希望每次显示步骤:科学模拟通常会保存每5000步左右.如果你想看漂亮的电影,你可以插入一些移动平均线以消除噪音(你的模拟没有温度或任何类型的东西,所以没有平均你可能会很好).此外,我不确定您的数据结构是否针对作业进行了优化,或者您是否遇到了缓存未命中问题.我不是一个真正的专家,所以也许其他人可能想要对此进行权衡.最后,考虑不使用pow,而是快速反平方根或类似方法.