串行代码
首先是简单的NBody代码,这里设置每个质点的质量相等。初始速度为0,初始位置均匀分布。
类似于这样 (16 * 16 大小 ):
这样可以抽象出代码循环逻辑为,这里N为总的质点数( 如上图,则为256 ),然后每一次update时,计算每个质点的受力(加速度)、速度、位置:
void update() {
for (int i = 0; i < N; i++)
force(i);
for (int i = 0; i < N; i++)
velocity(i);
for (int i = 0; i < N; i++)
position(i);
}
那么我们需要一些数据结构来保存数据,分别保存质点的速度、加速度、位置(2维的):
double* vBallvx;
double* vBallvy;
double* vBallax;
double* vBallay;
double* vBallpx;
double* vBallpy;
double r = 0.01; // 质点半径
const double GM = 6.67E-7; // G*M 用于计算加速度
const double delta_t = 0.01; // 表示每次update时经过的时间
void init(){
vBallvx = new double[N];
vBallvy = new double[N];
vBallax = new double[N];
vBallay = new double[N];
vBallpx = new double[N];
vBallpy = new double[N];
static int size = sqrt(N);
for (int i = 0; i < N; i++)
{
vBallpx[i] = 0.01 * (i % size);
vBallpy[i] = 0.01 * (i / size);
vBallvx[i] = 0;
vBallvy[i] = 0;
vBallax[i] = 0;
vBallay[i] = 0;
}
}
这样一个简单的串行多体问题代码就可以这样:
#include <iostream>
using namespace std;
double* vBallvx;
double* vBallvy;
double* vBallax;
double* vBallay;
double* vBallpx;
double* vBallpy;
int N = 256; // 这里N不设置为常数,并且vBall等数据动态分配是为了方便之后改为MPI
const double r = 0.01; // 质点半径
const double GM = 6.67E-7; // G*M 用于计算加速度
const double delta_t = 0.01; // 表示每次update时经过的时间
void init() {
vBallvx = new double[N];
vBallvy = new double[N];
vBallax = new double[N];
vBallay = new double[N];
vBallpx = new double[N];
vBallpy = new double[N];
static int size = sqrt(N);
for (int i = 0; i < N; i++)
{
vBallpx[i] = 0.01 * (i % size);
vBallpy[i] = 0.01 * (i / size);
vBallvx[i] = 0;
vBallvy[i] = 0;
vBallax[i] = 0;
vBallay[i] = 0;
}
}
void force(int index) {
vBallax[index] = 0;
vBallay[index] = 0;
for (int i = 0; i < N; i++) {
if (i == index) continue;
double dx = vBallpx[index] - vBallpx[i];
double dy = vBallpy[index] - vBallpy[i];
double d = dx * dx + dy * dy;
if (d < r * r) d = r * r; // ^(3/2)
d *= sqrt(d);
vBallax[index] += GM * (dx) / d;
vBallay[index] += GM * (dy) / d;
}
}
void velocity(int index) {
vBallvx[index] += vBallax[index] * delta_t;
vBallvy[index] += vBallay[index] * delta_t;
}
void position(int index) {
static int size = sqrt(N) - 1;
vBallpx[index] += vBallvx[index] * delta_t;
if (vBallpx[index] > size / 100.0) vBallpx[index] = size / 100.0;
if (vBallpx[index] < 0) vBallpx[index] = 0;
vBallpy[index] += vBallvy[index] * delta_t;
if (vBallpy[index] > size / 100.0) vBallpy[index] = size / 100.0;
if (vBallpy[index] < 0) vBallpy[index] = 0;
}
void update() {
for (int i = 0; i < N; i++)
force(i);
for (int i = 0; i < N; i++)
velocity(i);
for (int i = 0; i < N; i++)
position(i);
}
int main() {
init();
int i = 0;
while (true)
{
update();
cout << "迭代" << ++i << "次" << endl;
}
return 0;
}
OpenGL 可视化的代码就是这样,(使用GLUT这个OpenGL库):
#include <iostream>
#include <gl/glut.h>
using namespace std;
double* vBallvx;
double* vBallvy;
double* vBallax;
double* vBallay;
double* vBallpx;
double* vBallpy;
int N = 256; // 这里N不设置为常数,并且vBall等数据动态分配是为了方便之后改为MPI
const double r = 0.01; // 质点半径
const double GM = 6.67E-7; // G*M 用于计算加速度
const double delta_t = 0.01; // 表示每次update时经过的时间
void init() {
vBallvx = new double[N];
vBallvy = new double[N];
vBallax = new double[N];
vBallay = new double[N];
vBallpx = new double[N];
vBallpy = new double[N];
static int size = sqrt(N);
for (int i = 0; i < N; i++)
{
vBallpx[i] = 0.01 * (i % size);
vBallpy[i] = 0.01 * (i / size);
vBallvx[i] = 0;
vBallvy[i] = 0;
vBallax[i] = 0;
vBallay[i] = 0;
}
}
void force(int index) {
vBallax[index] = 0;
vBallay[index] = 0;
for (int i = 0; i < N; i++) {
if (i == index) continue;
double dx = vBallpx[index] - vBallpx[i];
double dy = vBallpy[index] - vBallpy[i];
double d = dx * dx + dy * dy;
if (d < r * r) d = r * r; // ^(3/2)
d *= sqrt(d);
vBallax[index] += GM * (dx) / d;
vBallay[index] += GM * (dy) / d;
}
}
void velocity(int index) {
vBallvx[index] += vBallax[index] * delta_t;
vBallvy[index] += vBallay[index] * delta_t;
}
void position(int index) {
static int size = sqrt(N) - 1;
vBallpx[index] += vBallvx[index] * delta_t;
if (vBallpx[index] > size / 100.0) vBallpx[index] = size / 100.0;
if (vBallpx[index] < 0) vBallpx[index] = 0;
vBallpy[index] += vBallvy[index] * delta_t;
if (vBallpy[index] > size / 100.0) vBallpy[index] = size / 100.0;
if (vBallpy[index] < 0) vBallpy[index] = 0;
}
void update() {
for (int i = 0; i < N; i++)
force(i);
for (int i = 0; i < N; i++)
velocity(i);
for (int i = 0; i < N; i++)
position(i);
}
// ----------------------------------------- OpenGL 函数分割线 ----------------------------------
void drawBall() {
glColor3f(0, 0, 0);
glPointSize(4);
glBegin(GL_POINTS);
for (int i = 0; i < N; i++)
glVertex2f(vBallpx[i] * 2, vBallpy[i] * 2);
glEnd();
}
void display(void) {
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
drawBall();
glutSwapBuffers();
}
void timeToEva(int ms) {
display();
glutTimerFunc(ms, timeToEva, ms);
}
void resize(int w, int h) {
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
glViewport(0, 0, w, h);
gluPerspective(45.0f, 1.0f * w / h, 1.0f, 100.0f);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
}
// 调用后阻塞
void initOpenGL(int argc, char* argv[], void (*update)()) {
glutInit(&argc, argv);
glutInitDisplayMode(GLUT_DEPTH | GLUT_DOUBLE | GLUT_RGBA);
glutInitWindowPosition(100, 100);
glutInitWindowSize(1000, 1000);
glutCreateWindow("Demo");
glutDisplayFunc(&display);
int ms = 30;
glutTimerFunc(ms, timeToEva, ms); // 每30ms执行一次,显示刷新(渲染)
glutIdleFunc(update); // CPU空闲时执行update
glClearColor(0.93f, 0.93f, 0.93f, 0.0f);
glEnable(GL_DEPTH_TEST);
glutMainLoop();
}
// ----------------------------------------- OpenGL 函数分割线 ----------------------------------
int main(int argc, char* argv[]) {
init();
initOpenGL(argc, argv, update);
return 0;
}
使用MPI代码
MPI有多种进程间通信函数,这里先选用的MPI_Bcast,然而一对多等集合操作没有tag,所以我们这里牺牲性能,用MPI_Send和MPI_Recv实现一个my_bcast。
参考: https://blog.csdn.net/a1066196847/article/details/89709720
参考: https://stackoverflow.com/questions/64124667/is-mpi-bcast-supposed-to-work-with-mpi-ibcast
代码如下:
void my_bcast(void* data, int count, MPI_Datatype datatype, int root, MPI_Comm communicator, int tag) {
int world_rank;
MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);
int world_size;
MPI_Comm_size(MPI_COMM_WORLD, &world_size);
if (world_rank == root) {
//如果是源进程的话,就要发送数据给到其他所有的进程
for (int i = 0; i < world_size; i++) {
if (i != world_rank) {
// 要发送count这么多的data,data的数据类型是datatype,要发送到i进程,数据的tag是0
MPI_Send(data, count, datatype, i, tag, communicator);
}
}
}
else {
//如果不是源进程,那么最终目的就只剩一个了,就是接收
MPI_Recv(data, count, datatype, root, tag, communicator, MPI_STATUS_IGNORE);
}
}
回到我们的多体问题来,时间主要耗在了计算加速度,即for循环那里,时间复杂度为O(N2),因此我们需要重点优化这里,思路就是将质点均匀的分配给每个进程。这里为了方便讨论,假设有4个进程,如下:
当计算加速度时,每个进程向其他所有进程广播当前管理的质点的位置。代码类似于这样,numprocs表示当前的进程数,如上图则为4,vBallpxAll、vBallpyAll表示同步的位置信息:
void update() {
for (int i = 0; i < numprocs; i++) {
my_bcast(vBallpxAll + i * N, N, MPI_DOUBLE, i, MPI_COMM_WORLD, 0);
my_bcast(vBallpyAll + i * N, N, MPI_DOUBLE, i, MPI_COMM_WORLD, 1); // TODO
}
for (int i = 0; i < N; i++)
force(i);
for (int i = 0; i < N; i++)
velocity(i);
for (int i = 0; i < N; i++)
position(i);
}
这里需要每个进程创建256大小的vBallpxAll、vBallpyAll,用来存放所有的质点的位置信息。而速度和加速度则不需要同步,只需要创建256/4 = 64大小的数组。每次迭代update前,将所有质点的位置信息进行同步。这样就可以正确的运行程序了。
最终代码如下:
#include <iostream>
#include <gl/glut.h>
#include<mpi.h>
using namespace std;
double* vBallvx;
double* vBallvy;
double* vBallax;
double* vBallay;
double* vBallpxAll;
double* vBallpyAll;
int N = 256; // 这里N不设置为常数,并且vBall等数据动态分配是为了方便之后改为MPI
const double r = 0.01; // 质点半径
const double GM = 6.67E-7; // G*M 用于计算加速度
const double delta_t = 0.01; // 表示每次update时经过的时间
int myid, numprocs; // 当前进程id 当前创建进程数
void my_bcast(void* data, int count, MPI_Datatype datatype, int root, MPI_Comm communicator, int tag) {
int world_rank;
MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);
int world_size;
MPI_Comm_size(MPI_COMM_WORLD, &world_size);
if (world_rank == root) {
//如果是源进程的话,就要发送数据给到其他所有的进程
for (int i = 0; i < world_size; i++) {
if (i != world_rank) {
// 要发送count这么多的data,data的数据类型是datatype,要发送到i进程,数据的tag是0
MPI_Send(data, count, datatype, i, tag, communicator);
}
}
}
else {
//如果不是源进程,那么最终目的就只剩一个了,就是接收
MPI_Recv(data, count, datatype, root, tag, communicator, MPI_STATUS_IGNORE);
}
}
void init() {
vBallvx = new double[N];
vBallvy = new double[N];
vBallax = new double[N];
vBallay = new double[N];
vBallpxAll = new double[N * numprocs];
vBallpyAll = new double[N * numprocs];
static int size = sqrt(N);
static int numprocsQ = sqrt(numprocs);
int offsetSize1 = myid % numprocsQ * size;
int offsetSize2 = myid / numprocsQ * size;
for (int i = 0; i < N; i++)
{
vBallpxAll[i + myid * N] = 0.01 * (i % size + offsetSize1);
vBallpyAll[i + myid * N] = 0.01 * (i / size + offsetSize2);
vBallvx[i] = 0;
vBallvy[i] = 0;
vBallax[i] = 0;
vBallay[i] = 0;
}
}
void force(int index) {
vBallax[index] = 0;
vBallay[index] = 0;
for (int i = 0; i < N * numprocs; i++) {
if (i == index + myid * N) continue;
double dx = vBallpxAll[index + myid * N] - vBallpxAll[i];
double dy = vBallpyAll[index + myid * N] - vBallpyAll[i];
double d = dx * dx + dy * dy;
if (d < r * r) d = r * r; // ^(3/2)
d *= sqrt(d);
vBallax[index] += GM * (dx) / d;
vBallay[index] += GM * (dy) / d;
}
}
void velocity(int index) {
vBallvx[index] += vBallax[index] * delta_t;
vBallvy[index] += vBallay[index] * delta_t;
}
void position(int index) {
static int size = sqrt(N * numprocs) - 1;
vBallpxAll[index + myid * N] += vBallvx[index] * delta_t;
if (vBallpxAll[index + myid * N] > size / 100.0) vBallpxAll[index + myid * N] = size / 100.0;
if (vBallpxAll[index + myid * N] < 0) vBallpxAll[index + myid * N] = 0;
vBallpyAll[index + myid * N] += vBallvy[index] * delta_t;
if (vBallpyAll[index + myid * N] > size / 100.0) vBallpyAll[index + myid * N] = size / 100.0;
if (vBallpyAll[index + myid * N] < 0) vBallpyAll[index + myid * N] = 0;
}
void update() {
for (int i = 0; i < numprocs; i++) {
my_bcast(vBallpxAll + i * N, N, MPI_DOUBLE, i, MPI_COMM_WORLD, 0);
my_bcast(vBallpyAll + i * N, N, MPI_DOUBLE, i, MPI_COMM_WORLD, 1); // TODO
}
for (int i = 0; i < N; i++)
force(i);
for (int i = 0; i < N; i++)
velocity(i);
for (int i = 0; i < N; i++)
position(i);
}
// ----------------------------------------- OpenGL 函数分割线 ----------------------------------
void drawBall() {
glColor3f(0, 0, 0);
glPointSize(4);
glBegin(GL_POINTS);
for (int i = 0; i < N * numprocs; i++)
glVertex2f(vBallpxAll[i] * 2, vBallpyAll[i] * 2);
glEnd();
}
void display(void) {
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
drawBall();
glutSwapBuffers();
}
void timeToEva(int ms) {
display();
glutTimerFunc(ms, timeToEva, ms);
}
void resize(int w, int h) {
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
glViewport(0, 0, w, h);
gluPerspective(45.0f, 1.0f * w / h, 1.0f, 100.0f);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
}
// 调用后阻塞
void initOpenGL(int argc, char* argv[], void (*update)()) {
glutInit(&argc, argv);
glutInitDisplayMode(GLUT_DEPTH | GLUT_DOUBLE | GLUT_RGBA);
glutInitWindowPosition(100, 100);
glutInitWindowSize(1000, 1000);
glutCreateWindow("Demo");
glutDisplayFunc(&display);
int ms = 30;
glutTimerFunc(ms, timeToEva, ms); // 每30ms执行一次,显示刷新(渲染)
glutIdleFunc(update); // CPU空闲时执行update
glClearColor(0.93f, 0.93f, 0.93f, 0.0f);
glEnable(GL_DEPTH_TEST);
glutMainLoop();
}
// ----------------------------------------- OpenGL 函数分割线 ----------------------------------
int main(int argc, char* argv[]) {
MPI_Init(&argc, &argv);
MPI_Comm_rank(MPI_COMM_WORLD, &myid); // pid
MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
N /= numprocs;
init();
if (myid != 0) {
while (true)
{
update();
}
}
else
initOpenGL(argc, argv, update);
return 0;
}
最后根据硬件属性,在相应的for循环前面添加#pragma omp parallel for
,让OpenMP自动进行共享内存的并行计算。整个多体问题的简单并行代码即完成。
附录 环境配置
OpenGL glut 配置
MPI 配置:
暂无