加快天文上多体模拟计算之探讨

1、加快天文上多體模擬計算之探討研究者：陳妙昕 薛淯勻指導者：闕志鴻 教授 曾耀寰 老師 李美英 老師壹、緒論一、研究動機參觀世貿電腦展覽時，發現大眾口中經常提起的超級電腦，並沒有如想像中那般的超級。記憶體是超級電腦的重要配備，限制了電腦執行程式運算的能力。所以我們希望能夠減少記憶體抓取不必要資料的次數，以達到節約資源，並做更有效的計算。二、研究目的繼Tree Code分格方法之後，為了使電腦在計算宇宙中兩兩質點間的引力大小時，不致浪費過多的記憶體在抓取資料上，我們試圖尋找一種最佳的方法，把範圍內的所有質點以一直線相連，使電腦在計算所有質點的數據時，能走一條最短的路徑。三、研究問題我們所要探討的問

2、題是使用何種方法可以加快天文上多體模擬的計算，而最終目的則是希望這一條路徑能將相近的質點快速地連結起來，而對於遠方的質點團，就利用其質心來計算其作用力。貳、文獻探討 多質點平行切割演算法（A Parallel Hashed Oct-Tree N-Body Algorithm）Michael S. Warreh & John K. Salmon一、 摘要作者最近設計並且實施在多質點間力的演算與NlogN成比例，其計算作用力的準確性受到限制，但可經由定義一個參數來做調整。作者使用一個切碎的表格來進行位置記憶的製圖，而不使用指示物來顯示位相。這樣的方法可讓多層次處理器更有效率的抓取資料。二、 介紹在

3、複雜的物理系統中，多質點的計算已成為一種基礎的工具，從一些基礎物理作用力（例如：重力、凡德瓦爾力）代表系統密度分類的多質點系統可遵循力學發展。在自然界中多質點計算為主要的統計數據（除了分子動力學演算被直接複製的物理系統。）更多質點可以提供更準確及更複雜的例子，因而能導出更準確或更複雜的結果。不幸的是，最小的準確度需要比現在計算資源所允許多很多的質點。在多質點計算中，因為交互作用發生在每一對質點之間，所以算數值得規模就漸近於N2，當能維持可接受的準確度時，許多的努力就花費在減少這樣計算的複雜性。一種處理方法是插入分解的格子，這項方法的缺點是力學相當的規模無法被模製。在三維空間中，這項對力學範圍的

4、限制，對許多的計算來說是不足的。另一個方法是將交互作用分成近和遠的部分。由於作用力的大小與粒子間的距離有關，因此當距離所造的力是可忽略的，這樣的力就可以被完全省略。然而，這樣有意義的錯誤很難被分析，當許多有意義的聚集系統中，將減低這個方法的效率。這一小部份的質點就會被放置在近的部分。多質點演算法的概念是當每一質點維持足夠的準確度時，將連續接近的質點逐步分割。經使用這種樹枝狀結構( tree )之後，一般說來這種方法陳述了階級組織方法的多質點系統。做出好的決定且刪除不精準的交互作用是此演算法成功與否的決定性因素。這樣的決定被一種多極接受標準（multipole acceptance criter

5、ion MAC）所控制。三、演算作者根據在d度空間所有點座標的對應關係，經轉換變成一組組數字( 以二進位的形式呈現 )，定義了一個key。這整個函數程式就包括了把點座標數字轉換成整數的過程，然後再將轉換好的d個整數插入一個key中。在這其中，作者並沒有對d的範圍加以限制，但實際上作者的研究主要還是設定在以三度空間為主，如此一來，由三個點座標經轉換所組成的key即可以一個64位元的整數或是一個32位元組的數完整地呈現出來。跳脫出以往所使用的平凡的函數系統，作者這次的函數是比較類似於Morton ordering。這種作業程式也是把系統中的每一個質點已一個key來代表。此外，作者還希望把這種key

6、的表現形式運用到樹枝狀結構( tree )上。為了要區別在樹枝狀結構中較高層以及較低層的質點，作者預先準備了一個額外的“1”給所有key中最重要的那個，這樣一來就可以呈現出在同一個key的範圍內所有較高層的質點，並且清楚的指出每一個質點的所在位置。一般來說，每一個key都會對應出一些複合的資料，這些data則多為描述此一晶格( cell )範圍內的相關物理資料，例如質心、質量等等。為了要使key和記憶體上儲存資料的地方能夠互相對應，作者就必須利用“切割表格( hash table )”。在執行這個切碎的系統，作者利用一個非常簡單的程式：AND的指令，並逐次找出最不明顯的數字( the leas

7、t significant bits)繼續進行分割。然而在run程式的過程中，電腦仍可能因為所被給予的指令不夠清晰，而使整個程式在執行上花費過多的時間及造成許多的錯誤，所以作者在防止錯誤的發生上是透過一個相連結鏈子( chaining )來解決的。將key的範圍劃分為有層次的結構以及找出最不明顯的數字的方法，經過作者的測試，的確能有效地減少錯誤( collision )的發生。當所有的key在其範圍內皆擁有最少經質點座標轉換後的數字，這個切割的程式就以近趨完美了。在tree中更高程度的節點可以多樣的方式來建構，最簡單的方式便是從核心開始裝填每一個粒子，然後越過部分架構的tree。當兩個粒子落在

8、同一格cell時格子就會再被更細的分割，使得tree的結構得以有更深入的階層，將這兩個粒子分開在不同的格子中。參、研究方法首先利用sin, cos的疊加來取亂數，代表某一空間中的質點。再定出質點的 x, y, z 軸座標。利用以下四種排序方法以決定路徑的走法：一、依各質點與原點的距離定序。二、依質點的 x, y, z 座標相乘後定序。三、將質點的 x, y, z 座標和相加。四、將質點的 x, y, z 座標分別開平方根後相加。 為了檢驗所選出的方法是否走出最佳的路徑，我們設定以一質點為中心，取其一定範圍內的質點：一、定義 w 為限定範圍內的總質點數。二、定義 s 為限定範圍內兩兩質點序號差的

9、和。三、比較各種排序方法執行出的 w 和 s 值。若 w 值愈大，s 值愈小，意味著近距離內的質點被定在相近的序號上，此路徑也就愈成功。四、利用 sum1（所有 w 的和）、sum2（所有 s 的和），可明顯比較出程式的優劣。肆、研究步驟一、撒質點：利用sinq, cosq 的疊加，可組成任意圖形。二、定座標：定出質點的 x, y, z 軸座標。三、路徑：以各種方法將質點串聯起來。四、排序：(一) 依各質點與原點的距離定序。 (二) 依質點的 x, y, z 座標相乘後定序。 (三) 將質點的 x, y, z 座標和相加。 (四) 將質點的 x, y, z 座標分別開平方根後相加。五、檢驗：(

10、一) 以一質點為中心，取其一定範圍內的質點。(二) 定義 w 為限定範圍內的總質點數。(三) 定義 s 為限定範圍內兩兩質點序號差的和。(四) 比較各種排序方法執行出的 w 和 s 值。若 w 值愈大，s 值愈小，意味著近距離內的質點被定在相近的序號上，此路徑也就愈成功。(五) 利用 sum1、sum2，可明顯比較出程式的優劣。 伍、研究設計 在仔細研讀過“ A Parallel Hashed Oct-Tree N-Body Algorithm ”論文以及指導老師的建議下，我們選擇利用C語言來撰寫程式，並在Linux的作業系統下著手進行所有的程式，而在繪圖方面則是運用Gnu plot的軟體來完

11、成繪圖的工作。以下是我們的最終程式：四組排序的比較#include #include #include #include #define PI 3.#define MAX 5000#define swap(a,b,t) (t)=(a); (a)=(b); (b)=(t);typedef struct float x; float y; float z; float d1; float d2; float d3; float d4; double t; int s1; int s2; int s3; int s4; int w1; int w2; int w3; int w4; value;vo

12、id sort1(value *arr, int left, int right) float pivot; int i,j; value temp; if(leftright) i=left; j=right+1; pivot=arrleft.d1; do do i+; while(arri.d1pivot); if(ij) swap(arri,arrj,temp); while(ij); swap(arrleft,arrj,temp); sort1(arr,left,j-1); sort1(arr,j+1,right); void sort2(value *arr, int left, i

13、nt right)float pivot; int i,j; value temp; if(leftright) i=left; j=right+1; pivot=arrleft.d2; do do i+; while(arri.d2pivot); if(ij) swap(arri,arrj,temp); while(ij); swap(arrleft,arrj,temp); sort2(arr,left,j-1); sort2(arr,j+1,right); void sort3(value *arr, int left, int right) float pivot; int i,j; v

14、alue temp; if(leftright) i=left; j=right+1; pivot=arrleft.d3; do do i+; while(arri.d3pivot); if(ij) swap(arri,arrj,temp); while(ij); swap(arrleft,arrj,temp); sort3(arr,left,j-1); sort3(arr,j+1,right); void sort4(value *arr, int left, int right) float pivot; int i,j; value temp;if(leftright) i=left;

15、j=right+1; pivot=arrleft.d4; do do i+; while(arri.d4pivot); if(ij) swap(arri,arrj,temp); while(ij); swap(arrleft,arrj,temp); sort4(arr,left,j-1); sort4(arr,j+1,right); int main(void) value vMAX; float a,b,c,d,e,f,t,m=1,n=1,o=1,p=1,q=1,r=1; int i,j,sum11,sum12,sum21,sum22,sum31,sum32,sum41,sum42; sra

16、nd(unsigned)time(NULL); freopen(test5.txt,w,stdout); for(i=0;iMAX;i+) a=(float)rand()/RAND_MAX)*0.5*PI; b=(float)rand()/RAND_MAX)*0.5*PI; c=(float)rand()/RAND_MAX)*0.5*PI; d=(float)rand()/RAND_MAX)*0.5*PI; e=(float)rand()/RAND_MAX)*0.5*PI; f=(float)rand()/RAND_MAX)*0.5*PI; vi.x=m*sin(a)+n*cos(b); vi

17、.y=o*sin(c)+p*cos(d); vi.z=q*sin(e)+r*cos(f); vi.d1=(vi.x)*(vi.x)+(vi.y)*(vi.y)+(vi.z)*(vi.z); vi.d2=(vi.x)*(vi.y)*(vi.z); vi.d3=(vi.x)+(vi.y)+(vi.z); vi.d4=sqrt(vi.x)+sqrt(vi.y)+sqrt(vi.z); / 排第一次 sort1(v,0,MAX-1); for(i=0;iMAX;i+) vi.s1 = 0; vi.w1 = 0; vi.t = 0.; for(j=0;jMAX;j+) if(i!=j) t=(vi.x-

18、vj.x)*(vi.x-vj.x)+(vi.y-vj.y)*(vi.y-vj.y)+(vi.z-vj.z)*(vi.z-vj.z); if(t0.1) vi.w1 = vi.w1 + 1; vi.s1 = vi.s1+fabs(i-j); vi.t += t; sum11 = 0; sum12 = 0; for(i=0;iMAX;i+) sum11+=vi.w1; sum12+=vi.s1; printf(sum11=%d sum12=%dn,sum11,sum12); for(i=0;iMAX;i+) printf(x=%f y=%f z=%f t=%f w1=%d s1=%d d1=%fn

19、,vi.x, vi.y, vi.z, vi.t, vi.w1, vi.s1, vi.d1); printf(nnn);/ 排第二次 sort2(v,0,MAX-1); for(i=0;iMAX;i+) vi.s2 = 0; vi.w2 = 0; vi.t = 0.; for(j=0;jMAX;j+) if(i!=j) t=(vi.x-vj.x)*(vi.x-vj.x)+(vi.y-vj.y)*(vi.y-vj.y)+(vi.z-vj.z)*(vi.z-vj.z); if(t0.1) vi.w2 = vi.w2 + 1; vi.s2 = vi.s2+fabs(i-j); vi.t += t; s

20、um21 = 0; sum22 = 0; for(i=0;iMAX;i+) sum21+=vi.w2; sum22+=vi.s2; printf(sum21=%d sum22=%dn,sum21,sum22); for(i=0;iMAX;i+) printf(x=%f y=%f z=%f t=%f w2=%d s2=%d d2=%fn,vi.x, vi.y, vi.z, vi.t, vi.w2, vi.s2, vi.d2); printf(nnn);/ 排第三次 sort3(v,0,MAX-1); for(i=0;iMAX;i+) vi.s3 = 0; vi.w3 = 0; vi.t = 0.

21、; for(j=0;jMAX;j+) if(i!=j) t=(vi.x-vj.x)*(vi.x-vj.x)+(vi.y-vj.y)*(vi.y-vj.y)+(vi.z-vj.z)*(vi.z-vj.z); if(t0.1) vi.w3 = vi.w3 + 1; vi.s3 = vi.s3+fabs(i-j); vi.t += t; sum31 = 0; sum32 = 0; for(i=0;iMAX;i+) sum31+=vi.w3; sum32+=vi.s3; printf(sum31=%d sum32=%dn,sum31,sum32);for(i=0;iMAX;i+) printf(x=%

22、f y=%f z=%f t=%f w3=%d s3=%d d3=%fn,vi.x, vi.y, vi.z, vi.t, vi.w3, vi.s3, vi.d3); printf(nnn); / 排第四次 sort4(v,0,MAX-1); for(i=0;iMAX;i+) vi.s4 = 0; vi.w4 = 0; vi.t = 0.; for(j=0;jMAX;j+) if(i!=j) t=(vi.x-vj.x)*(vi.x-vj.x)+(vi.y-vj.y)*(vi.y-vj.y)+(vi.z-vj.z)*(vi.z-vj.z); if(t0.1) vi.w4 = vi.w4 + 1; v

23、i.s4 = vi.s4+fabs(i-j); vi.t += t; sum41 = 0; sum42 = 0; for(i=0;iMAX;i+) sum41+=vi.w4; sum42+=vi.s4; printf(sum41=%d sum42=%dn,sum41,sum42); for(i=0;iMAX;i+) printf(x=%f y=%f z=%f t=%f w4=%d s4=%d d4=%fn,vi.x, vi.y, vi.z, vi.t, vi.w4, vi.s4, vi.d4); 四種排序方法比較(取500個點)陸、結論 關於” Traveling Salesman Probl

24、em ” & “ nearest neighborhood method ”已在數學界存在悠久，且一直是困擾眾多數學家而至今仍無解的問題。我們在此僅提出四種簡單方法，爾後還能利用其他路徑來比較各路徑的優劣，而且不需要大量的知識背景。此外，在我們的實驗結論之中，利用相乘的方法是最佳的路徑，關於其背後所蘊含的數學理論，也是值得去探討、研究的，但這個部分就有賴於比較多的專業知識。雖然在上面的四個圖中，由於太多的線條使得我們難以分辨程式的優缺，但經由檢驗程式的計算，就能夠簡單的用數值來判斷。且在多次試驗之後，都由相承的方法奪得后冠，而開平方根後相加則是眾所期待的榜眼，但這一切仍有待嚴謹的證明。柒、討論與建議 經過一年來的鍛鍊，我們深深覺得看事情真的不能只看表面。當初選擇這個實驗室完全沒考慮到天文物理在實際實驗操作上的困