【发布时间】：2015-02-13 07:34:52
【问题描述】：
我正在尝试实现二维 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;
}
}
}
【问题讨论】：
-
您的
j循环似乎正在访问尚未由您的i循环初始化的数组元素。也许在使用i和j进行双循环之前运行整个i循环进行初始化？ -
降低代码的复杂度，减少出错的可能性。
-
@user3558391：我怀疑这里的问题至少部分是由于您使用的集成方法不准确（一阶Euler method）。您是否考虑过改用Verlet integration？当粒子相互坍缩时，它会在势能中达到奇点，因此那里可能会出现较大的数值误差。
-
我知道，这就是我的意思。在担心 OpenMP 之前，先让顺序版本正常工作。
-
我制作了以下输入文件：
2 100000000 500 500 0 0 0 0 1 530 500 0 0.005 0 0它的行为与我期望的一样，其中一颗行星围绕另一颗行星运行。现在，如果0.005被0取代，那么较轻的行星将与较大的行星相撞。这就是说，力最终将变得几乎无限，导致爆炸（参见例如gromacs.org/Documentation/Terminology/Blowing_Up）。在短距离内添加排斥力（比如 Lennard-Jones 类型），事情可能会看起来更好。
标签： c++ openmp simulation physics numerical-methods