如何實現在計算機上模擬宇宙的演變過程？

11-28

建立模型，輸入必要的參數，還有合適的計算機語言？
-----
本題已收錄至知乎圓桌 ?宇宙那麼大，更多討論歡迎關注。

應要求已加圖！！！（小心流量~）這個回答較長，另歡迎摳錯

首先回答題主的問題：
模型：最簡單的是基本的牛頓引力問題，即解泊松方程：

?^2 φ(x) = 4πGρ(x)

必要的參數：有太多，比如你要模擬的宇宙尺度的大小，粒子（就是天體）數目，等等。。。還有一些計算機的參數，通常要把這些參數全部寫進一個或幾個文件才可以。
計算機語言：C/C++，Fortran

然後展開了說：
這類問題其實就是多體問題，給出粒子的初始位置與速度，算出粒子與粒子之間的牛頓引力，計算其根據時間的演化。只是在此過程中，你需要考慮宇宙膨脹的效應。
當然，當粒子數量很少的時候，我們還是很容易計算的。對於計算機來說，也沒有任何負擔，但是當例子數目增加到成百上千萬的時候，我堅信你自己的計算機肯定是沒戲了~
那好，我們這樣做一個計算：
我們假設每個粒子質量一樣，只計算牛頓引力，因此每一個粒子的性質可以用它的位置(x1,x2,x3)和速度(v1,v2,v3)來描述，我們做如下假設：
每個粒子用單精度存儲，則：
2-Body問題需要佔用4B*6*2=48B的內存，計算1個粒子-粒子作用力（pp-force），可以解出Kepler 軌道。3-Body問題需要佔用4B*6*3=72B的內存，計算3個粒子-粒子作用力（pp-force），軌道沒有解析解。類推，N-body問題需要佔用4B*6*N=(24*N)B，計算(N-1)N/2個粒子-粒子作用力（pp-force），也沒有軌道的解析解。
這裡解釋一下，我們在這裡假設所有粒子分布在一個正方形「盒子」中，這個盒子滿足周期性邊界條件，也就是說粒子從一個邊出去，就從對應的另一邊進來。上圖 np是指初始我們撒粒子時，每一條邊撒的數目，於是np^3就是粒子總數了。我們把粒子的位置和速度信息存入xv.dat文件，可以看看這個文件的大小。假設你的計算機內存是4G，你估算一下，還是假設粒子信息以4B的精度存儲，你只能算約500^3的粒子數目。因此你需要大型計算機。我們繼續估算，當粒子數目很多的時候，如10^12兩級，無論什麼樣子的計算機也承受不了10^24這個數量級的pp-force的計算。（給一個圖來感受一下）

因此我們需要簡化這個計算過程，通常來說，有兩種方法：
1、Tree-Method（樹形法）：
也就是根據粒子密度不同，劃分網格，將粒子歸算與一個網格之內，計算格點格點與格點間的作用力。著名的Gadget程序就是這個方法，稍後在多科普一些Gadget，先回到pp-force的解決辦法上。

2、Particle-Mesh：
這種方法不同於樹形法，這個方法的網格大小是相同的，固定的，不會根據粒子密度不同而改變。這種方法下面還要分為 Cloud-in-Cell （CIC）和 Nearest-Grid-Particle （NGP）兩種解Poisson方程的演算法。
格點之間之間，計算格點與格點的作用力，格點內部，計算粒子與粒子的作用力，通過這種方法減少計算量。
在CubeP3M這個程序中，粗格點計算使用CIC，細格點使用NGP

如上圖，左面是粗格點，右面是細格點。
給一個CubeP3M的鏈接：
CubePM - CITA Computing
這個程序我一會也會說一些。

其實這裡面有一個小問題，就是假設粒子與粒子距離過於緊密，那麼根據平方法比例，粒子與粒子作用力就會很大很大，粒子加速度會很大很大很大很大，最後粒子就飛了。這會出問題的。這會兒就需要一個 force softening，即力的軟化問題。即在一個粒子的一定範圍內，不計算粒子與粒子之間的作用力。但是，即使簡化了這些計算，我們的小PC機還是不能滿足計算需求（除非你費大量時間一個一個地計算作用力），因此我們需要程序進行並行計算：Message Passing Interface（MPI），Open Multi-Processing （OpenMP），等等。。。

隨後就可以說一下計算的流程：
第一：生成初始條件（物質密度分布），即生成粒子的在「盒子」的位置
粒子的初始位置是根據初始功率譜確定的，可以通過CAMB計算轉移方程得到，而且我們計算得到的這個密度場是一個隨機高斯場。在操作時，我們先把粒子均勻擺在格點內，然後根據這個密度場將粒子進行移動，到它的初始位置。
第二：生成了計算每一步的作用力，然後算出粒子接下來的位置和速度，以此循環，知道模擬至當今宇宙（紅移等於0）。這一步說起來很簡單，但是實際上是程序中最複雜的一部分，核心是對泊松方程進行傅里葉變換，在傅里葉空間求解後返回實空間。（具體的看文獻吧，較複雜，太長不贅述，看Gadget和CubeP3M的文獻即可）當然，如果我們在這一步裡面不只有引力作用，比如我們加了磁流體加了重子，我們便需要考慮引力之外的物理過程，這會使得計算更複雜。
第三：計算物質密度功率譜（最主要），與觀測對比，檢驗結果。
這裡實際上就是計算物質密度分布的統計信息，根據尺度得到物質密度的功率譜。

（功率譜大概就長上圖這樣，實線是觀測，虛線是幾次模擬的）
算到這裡基本就結束了以此模擬運算，然後會有數據處理，不贅述。

除了純暗物質模擬，我們還可以加入重子（baryonic particles）過程，加入磁流體過程，等等。

下面舉例開始：
這中模擬宇宙演化的程序其實還是不少的。
通常來講，我們最常用的還是Gadget，而且很多程序都是基於它做的：Cosmological simulations with GADGET 文獻：GADGET: a code for collisionless and gasdynamical cosmological simulations，The cosmological simulation code GADGET-2，這個程序你在你自己的筆記本上也能玩的起來，但是規模很小而已~（下圖為當時自己筆記本跑的，因為筆記本壞了一次，參數都丟了。。。）

還有CubeP3M，這是我接觸的第一個這種程序http://arxiv.org/abs/1208.5098，jharno/cubep3m · GitHub這個程序目前的記錄是3萬億粒子（暗物質+中微子，應該是粒子數目的最高紀錄），在天河2號上運行的。目前數據正在處理中（下圖只是一個暗物質+中微子的小規測試時的模模擬）。

當然還有很多組在做這類工作，比如去年很火的Illustris：Properties of galaxies reproduced by a hydrodynamic simulation。如下圖（很帥！）：

國內也有做的，比如盤古模擬，

見下面這篇文獻：現代宇宙學中的數值模擬技術和應用--《中國科學:物理學力學天文學》2013年06期（這個講的應該很詳細了，中文的，把模擬原理很學術化地講一下，講的也很清楚）
最後做一個對未來的展望，數值模擬會怎樣發展呢？引用文獻的圖一給了很好的答案。粒子的數量級和年份呈線性增長。

最後，LZ確實姿勢水平有限，寫的比較淺顯，先這樣吧（也許還會有更新）

歡迎各位大神給補充，我也及時完善答案。另歡迎拍磚挑錯。

http://www.mpa-garching.mpg.de/gadget/

推薦universe sandbox引力模擬軟體

建立最小單元的特徵模型和行為模型，以及散布規律和初始狀態，然後就可以按下「宇宙運行」的開始按鈕啦！
下面是我寫的制定了物質密度，分布密度，分布範圍，初始速度的情況下，遵循萬有引力定律以及動量守恆定律（其實碰撞是由物理引擎實現的，這裡看作是遵循動量守恆定律）進行的宇宙物質模擬。發光點表示碰撞。

最後經過一定時間，形成了比較穩定的聚合體，可以認為是恆星，周邊圍繞核心做橢圓運動的可以看作行星，於是一個太陽系就誕生啦！

模擬過程中的感想：1、物質初速不能過快，否則靠引力無法聚合。2、物質密度不能過大，否則在聚合運動過程中會導致運動速度過快（畢竟是引力下的加速運動），仍然無法聚合。3、能形成穩定的太陽系時，參數範圍是比較小的，所以我們今天的太陽系能成為這個樣子，其實是蠻偶然的。4、我沒能模擬出黃道平面，行星軌道面成隨機分布，這個還請大牛指教一二。
5、其實這個模型稍微改改就是帶電粒子子群行為模型，去掉電荷就是分子碰撞模型，只是具體參數不同罷了，所以宏觀和微觀在很多方面真是很相近。

劉慈欣的《鏡子》是個思路而已在整個宇宙最簡單的時候開始模擬也就是大爆炸奇點That all started with the big bang!

你們設備爛爆了。

只要我們觀測地球和整個宇宙的現象進行反推，就能夠得到一個結論，我們的宇宙是被一種極其巨大的神秘能量模擬和有目的演進的結果。

就憑地球這點資源恐怕不夠啊