c++ - C++中的N体模拟

Question

我正在尝试实现二维 n-body 模拟的 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) 或其他位置，其上的力为零，这意味着它永远不会返回并且现在也具有相当大的速度。您的模拟有 blown up .

正如您所说，这与牛顿定律不符，因此这表明您的数值方案失败了。不用担心，有几种方法可以解决此问题。您已经预料到了，就像您添加了 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 interaction包含由泡利不相容原理引起的排斥。请注意，您可以增加排斥的锐度(如果您愿意，可以将 2.5 更改为 12.5)，这意味着排斥对远处的影响较小(好)，但它需要更准确地解决，导致 period 更小(不好)。理想情况下，您将从不会导致冲突的初始配置开始，但这几乎是不可能预测的。最后，您可能希望使用上面列出的方法的组合。

使用 OpenMP 可能会稍微加快速度，但您确实应该使用用于远程力的算法，例如 Barnes-Hut simulation .更多信息，请参见 2013 Summer School on Fast Methods for Long Range Interactions in Complex Particle Systems和 booklet (免费提供)他们发布了最新进展。您也可能不想显示每个时间步:科学模拟通常每 5000 步左右保存一次。如果你想要好看的电影，你可以插入一些移动平均线来消除噪音(你的模拟没有温度或类似的东西，所以你可能不用平均就可以了)。此外，我不确定您的数据结构是否针对该作业进行了优化，或者您是否可能遇到缓存未命中。我不是真正的专家，所以也许其他人可能想对此进行权衡。最后，考虑不要使用 pow，而是使用 fast inverse square root或类似的方法。