利用MPI解决N体问题

参考 MPI并行计算模拟N体问题

实验题目

题目：
利用MPI解决N体问题
实验内容：
N体问题是指找出已知初始位置、速度和质量的多个物体在经典力学情况下的后续运动。在本次实验中，你需要模拟N个物体在二维空间中的运动情况。通过计算每两个物体之间的相互作用力，可以确定下一个时间周期内的物体位置。
在本次实验中，初始情况下，N个小球等间隔分布在一个正方形的二维空间中，小球在运动时没有范围限制。每个小球间会且只会受到其他小球的引力作用。小球可以看成质点。小球移动不会受到其他小球的影响(即不会发生碰撞，挡住等情况)。你需要计算模拟一定时间后小球的分布情况，并通过MPI并行化计算过程。
实验要求
1. 有关参数要求如下：
  - 引力常数数值取6.67*10^11
  - 小球重量都为10000kg
  - 小球间的初始间隔为1cm，例：N=36时，初始的正方形区域为5cm*5cm
  - 小球初速为0
  - 其他未定义的参数均可自行拟定
2. 你的程序中，应当实现下面三个函数：
  - compute force(): 计算每个小球受到的作用力
  - compute velocities(): 计算每个小球的速度
  - compute positions(): 计算每个小球的位置
  典型的程序中，这三个函数应该是依次调用的关系。如果你的方法中不实现这三个函数，应当在报告中明确说明，并解释你的方法为什么不需要上述函数的实现。
3. 报告中需要有 N=64 和 N=256 的情况下通过调整并行度计算的程序执行时间和加速比

实验环境

操作系统	编译器	硬件配置
Ubuntu 16.04	mpicc	双核 4G内存

算法设计与分析

初始化条件

初始情况下，N个小球等间隔分布在一个正方形的二维空间中，小球在运动时没有范围限制。每个小球间会且只会受到其他小球的引力作用。小球可以看成质点。小球间的初始间隔为1cm

小球结构体

px和py为位置
ax和ay为加速度的水平和竖直分量
vx和vy为速度的水平和竖直分量

typedef struct object
{double px, py;double ax, ay;double vx, vy;
}object;

计算作用力(加速度)

d为两球间距离。R为小球半径，d的最小值为2R。

void compute_force(object *list, int size, int i)
{list[i].ax = 0;list[i].ay = 0;for(int k = 0; k < size; k++){if(k == i) continue;double dx = list[k].px - list[i].px;double dy = list[k].py - list[i].py;double d = sqrt(dx*dx + dy*dy);if(d < 2*R) d = 2*R;list[i].ax += G*M*dx/pow(d, 3);list[i].ay += G*M*dy/pow(d, 3);}
}

计算速度

delta_t是时间间隔，设置为0.001。理论上，该值越小，计算越精确。

void compute_velocity(object *list, int i)
{list[i].vx += list[i].ax * delta_t;list[i].vy += list[i].ay * delta_t;
}

计算位置

delta_t是时间间隔，设置为0.001。理论上，该值越小，计算越精确。

void compute_position(object *list, int i)
{list[i].px += list[i].vx * delta_t;list[i].py += list[i].vy * delta_t;
}

并行设计

初始化MPI环境，获取并行环境参数(总线程数、本地进程编号等)。

MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
MPI_Comm_rank(MPI_COMM_WORLD, &myid);

从0号进程获取输入数量和周期，并将其广播出去。

if(myid == 0){printf("数量:");scanf("%d", &n);l = (int)sqrt((double)n);printf("周期:");scanf("%d", &t);start = MPI_Wtime();
}
MPI_Bcast(&n, 1, MPI_INT, 0, MPI_COMM_WORLD);   //将n值广播出去
MPI_Bcast(&t, 1, MPI_INT, 0, MPI_COMM_WORLD);   //将t值广播出去
MPI_Bcast(&l, 1, MPI_INT, 0, MPI_COMM_WORLD);   //将l值广播出去

初始化列表。

object *list = (object*)malloc(n * sizeof(object));
for(int i = 0; i < n; i++){list[i].vx = 0;list[i].vy = 0;list[i].px = i%l*0.01;list[i].py = i/l*0.01;
}

核心循环

由于计算加速度时需要使用位置信息，所以更新加速度时不能同时更新位置信息，需要放在不同的循环中。
每个进程需要进程间通信以获取计算所需的位置信息。

for(int j = 0; j < t; j++){for(int i = n/numprocs*myid; i < n/numprocs*(myid+1); i++){compute_force(list, n, i);}MPI_Barrier(MPI_COMM_WORLD);for(int i = n/numprocs*myid; i < n/numprocs*(myid+1); i++){compute_velocity(list, i);compute_position(list, i);}MPI_Barrier(MPI_COMM_WORLD);}

结果统计

N=64 初始状态

N=64 1000周期

N=64 2000周期

N=256 初始状态

N=256 1000周期

N=256 2000周期

规模/线程数/运行时间(s)	1	2
N=64 1000周期	0.807149	0.492253
N=64 2000周期	1.612149	0.964136
N=256 1000周期	12.473251	7.423541
N=256 2000周期	24.889742	14.690790

规模/线程数/加速比	1	2
N=64 1000周期	1	1.640
N=64 2000周期	1	1.672
N=256 1000周期	1	1.680
N=256 2000周期	1	1.694

分析与总结

实验中没有4线程和8线程的数据，主要是因为实验时发现这两种情况下运行时间非常久，严重不符合预期。至于原因，暂未找到。
时间间隔delta_t表示每个周期的间隔，其值会影响小球分布形状。当delta_t=0.001时，小球分布图如上图所示，当delta_t=0.01时，小球分布图更趋于圆形。

利用MPI解决N体问题相关推荐

利用交换机解决局域网ARP问题(51CTO博客出书活动)
此文参加"51CTO博客出书:IT技术案例大征集"活动,欢迎讨论,谢绝转载! *案例名称: <利用交换机解决局域网ARP问题> *技术范围: 交换 ...
利用多线程解决多业务不同定时区间歇触发问题的一种方法
利用多线程解决多业务不同定时区间歇触发问题的一种方法参考文章: (1)利用多线程解决多业务不同定时区间歇触发问题的一种方法 (2)https://www.cnblogs.com/naaoveGIS/ ...
利用ComponentWillReceiveProps解决异步问题
利用ComponentWillReceiveProps解决异步问题参考文章: (1)利用ComponentWillReceiveProps解决异步问题 (2)https://www.cnblogs. ...
利用axios解决跨域的问题
利用axios解决跨域的问题参考文章: (1)利用axios解决跨域的问题 (2)https://www.cnblogs.com/Irelia/p/9972990.html 备忘一下.
深入跨域问题(2) - 利用 CORS 解决跨域
阅读目录: 深入跨域问题(1) - 初识 CORS 跨域资源共享: 深入跨域问题(2) - 利用 CORS 解决跨域(本篇) 深入跨域问题(3) - 利用 JSONP 解决跨域深入跨域问题(4) - ...
2.3 利用正规化解决过拟合问题-机器学习笔记-斯坦福吴恩达教授
利用正规化解决过拟合问题在之前的文章中,我们认识了过拟合问题,通常,我们有如下策略来解决过拟合问题: 减少特征数,显然这只是权宜之计,因为特征意味着信息,放弃特征也就等同于丢弃信息,要知道,特征的获 ...
MAT之ACA：利用ACA解决TSP优化最佳路径问题
MAT之ACA:利用ACA解决TSP优化最佳路径问题目录输出结果实现代码输出结果实现代码 load citys_data.mat n = size(citys,1); D = zeros(n ...
利用 Arthas 解决启动 HDFS StandbyNameNode 加载 EditLog 慢的问题
利用 Arthas 解决启动 StandbyNameNode 加载 EditLog 慢的问题公司新搭 HDFS 集群,namenode做ha,但是在启动 StandbyNamenode 节点的时候出 ...
c++获取输入数字的位数/获取位数并且将其存入数组中/获取位数存入数组并且利用它解决实际问题
1.仅仅获取数字的位数: int checkNumber(int num){int weishu = 0;for (int i = 1, num1 = num, shuzu = 0; num / i; ...