如何用python模擬一個星系？

11-27

用python這種不擅長硬體計算的語言真不合適。一般都是用C，最富盛名的模擬軟體是Gadget-2. 不過一般這個軟體是跑在伺服器上的，個人pc如果裝4路E5倒也應該跑的不錯。

如果你有興趣進行模擬，可以參照以下網址:
https://astrobites.org/2011/04/02/installing-and-running-gadget-2/

效果可參見常峻的視頻，其實這個視頻可以部分說明為什麼銀河系和仙女星系碰撞之後會有那麼多奇妙的天空背景。

如果某一天 M31 和銀河系相互併合，在地球的夜晚長期觀測會是一副怎樣的景象？

在一開始,還是要感謝一下凌晨曉驥.

我的這個回答相當於對他的答案的一個補充.主要目標是幫助大家安裝Gadget-2,並引導你們自己畫出銀河碰撞的視頻.

那麼開始了!

首先看一下我們的目標

星系碰撞模擬_野生技術協會_科技_bilibili_嗶哩嗶哩bilibili.com視頻

我用的是一台AMD ryzen1600的台式機,安裝了windows 1709系統.我會將Gadget-2安裝在wsl中.至於怎麼開啟wsl的內容我就不多介紹了,請自行百度一下.注意這個wsl還要有x11環境,windows裡面也要安裝x11伺服器,推薦使用VcXsrv.

我們假設gcc vim zlib wget等等常用的依賴我們都安裝好了現在只差安裝與Gadget-2緊密相關的MPI庫 GSL庫 FFTW庫和HDF5庫就可以安裝Gadget-2了,缺少其他庫可以根據報錯自行安裝依賴.

安裝MPI庫

Gadget-2需要MPI提供並行計算的支持,而MPI有幾種不同的實現.這裡我們選擇安裝MPICH,它是MPI庫的一種主流實現.因為我們的環境是wsl,也就是ubuntu.所以可以很方便的運行下面的命令來安裝MPICH.

sudo apt install mpich

安裝GSL庫

GSL全稱GNU scientific library,我們通過下面的命令安裝.

wget https://mirror.tuna.tsinghua.edu.cn/gnu/gsl/gsl-1.9.tar.gz tar -xzf gsl-1.9.tar.gz cd gsl-1.9 ./configure --prefix=/usr/local make sudo make install cd

安裝FFTW庫

FFTW是一個快速傅立葉變換庫,最新版本是FFTW3系列,但由於FFTW3去掉了對MPI的支持,所以我們只能安裝FFTW2的最後一個版本以保持兼容性.

wget http://www.fftw.org/fftw-2.1.5.tar.gz tar -xzf fftw-2.1.5.tar.gz cd fftw-2.1.5 ./configure --enable-mpi --enable-type-prefix --enable-flaot --prefix=/usr/local make sudo make install cd

安裝HDF5

如果你不想自己寫Gadget的文件格式parser的話還是把HDF5庫安裝上吧.由於兼容性問題,我們就不安裝最新版本的HDF5了.

wget https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.6/hdf5-1.6.9/src/hdf5-1.6.9.tar.gz tar -xzf hdf5-1.6.9.tar.gz cd hdf5-1.6.9 ./configure --prefix=/usr/local make sudo make install cd

安裝Gadget

前面的準備工作完成了,我們現在開始安裝Gadget.

wget http://wwwmpa.mpa-garching.mpg.de/gadget/gadget-2.0.7.tar.gz tar -xzf gadget-2.0.7.tar.gz cd Gadget-2.0.7/Gadget2 ls

這時我們列出了Gadget2文件夾下面所有的文件,然後用任何編輯器編輯其中的Makefile文件.文件的內容大致如下:

首先我們只需要確保下面這幾行被注釋掉就可以了.

#OPT += -DPERIODIC #OPT += -DPMGRID=128

然後我們要修改第101行左右的一部分內容.從

ifeq ($(SYSTYPE),"MPA") CC = mpicc OPTIMIZE = -O3 -Wall GSL_INCL = -I/usr/common/pdsoft/include GSL_LIBS = -L/usr/common/pdsoft/lib -Wl,"-R /usr/common/pdsoft/lib" FFTW_INCL= FFTW_LIBS= MPICHLIB = HDF5INCL = HDF5LIB = -lhdf5 -lz endif

修改為

ifeq ($(SYSTYPE),"MPA") CC = mpicc OPTIMIZE = -O3 -Wall GSL_INCL = -I/usr/local/include GSL_LIBS = -L/usr/local/lib -Wl,"-R /usr/local/lib" FFTW_INCL= -I/usr/local/include FFTW_LIBS= -L/usr/local/lib MPICHLIB = -L/usr/local/lib HDF5INCL = -I/usr/local/include HDF5LIB = -L/usr/local/lib -lhdf5 -lz endif

保存文件然後執行下面的命令

make

現在,Gadget-2已經成功安裝了.

下面進入第二個部分,測試模擬過程.

執行下面的命令.

cd cd Gadget-2.0.7 mkdir galaxy cp Gadget2/Gadget2 galaxy/ cp Gadget2/parameterfiles/galaxy.param galaxy/ cp ICs/galaxy_littleendian.dat galaxy/ cd galaxy mkdir snaps

然後編輯剛剛拷貝過來的galaxy.param文件,改動如下:

將前幾行的

% Relevant files

InitCondFile ICs/galaxy_littleendian.dat OutputDir galaxy/

改為

% Relevant files

InitCondFile ./galaxy_littleendian.dat OutputDir ./snaps

保存並執行

mpirun -np 4 ./Gadget2 galaxy.param

如果屏幕不斷快速閃過一堆字元,說明你成功了,等待結束即可.

但由於我個人並不是太喜歡原有的輸出格式,所以我們還可以繼續修改一下galaxy.param文件

找到第27行SnapFormat,將值修改成3.

第37行TimeMax,將值修改成10.

第47行TimeBetSnapshot,將值修改成0.01.

然後運行

mpirun -np 4 ./Gadget2 galaxy.param

等待運行結束,我們就在snap文件夾得到了大量的snapshot文件.

和之前的不一樣這次的snapshot文件的拓展名是hdf5.

現在我們獲得了數據文件,離把圖畫出來就差一步了.

關於畫圖,我們需要安裝python python-tk 和一些其他的python庫,你要是覺得麻煩找一台普通的linux或者mac也行.

sudo apt install python python-tk python-pip sudo pip isntall numpy matplotlib tables

然後在galaxy文件夾下面創建兩個文件和一個文件夾

touch analysis.py touch job.sh chmod a+x job.sh chomd a+x analysis.py mkdir png

然後把下面這兩個腳本的內容分別拷貝進analysis.py和job.sh裡面,腳本寫的太爛不要笑...

第一個是analysis.py

#!/usr/bin/env python # -*- coding:utf-8 -*-


from sys import argv

import numpy as np

import tables

import time

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D
NAME = argv[1]

I = int(argv[2])

OUTPUT = argv[3]
fig = plt.figure(figsize=(16, 12), dpi=80, facecolor=(0.0, 0.0, 0.0))

ax = Axes3D(fig)

h5 = tables.open_file(NAME % I)

halo = h5.root.PartType1

disk = h5.root.PartType2

ax.grid(False)

ax.set_axis_off()

ax.set_xlim(-200, 200)

ax.set_ylim(-200, 200)

ax.set_zlim(-200, 200)
ax.scatter(disk.Coordinates[:,0],

           disk.Coordinates[:,1],

           disk.Coordinates[:,2],

           color="black", s=0.5)

ax.scatter(halo.Coordinates[:,0],

           halo.Coordinates[:,1],

           halo.Coordinates[:,2],

           color="blue", s=0.01)

plt.savefig("%s%d" % (OUTPUT, I))

h5.close()

plt.close("all")

print "第%d個圖片完成了" % I

第二個是job.sh

#!/usr/bin/env bash


if [ $# -ne 4 ]

then

        echo "用法:"

        echo "      ./job.sh ./analysis.py snapshot的路徑 snapshot的數量 保存圖片的文件夾"

        echo "例子:"

        echo "      ./job.sh ./analysis.py ./snap/snapshot_%03d.hdf5 3000 ./png/snap_"

        echo "輸入文件中%03d是替換文件名包含的數字的,輸出文件名只是一半,腳本會自動在後面加上數字和拓展名"

fi
analysisScript=$1

inputPath=$2

cycleIndex=$3

outputPath=$4

echo "正在執行請稍候..."

maxProcess=10

fifoFile="/tmp/$$.fifo"
mkfifo $fifoFile

exec 6&<&>$fifoFile rm -rf $fifoFile for((i=1; i&<=$maxProcess; i++)); do echo done&>6 for index in `seq 0 $cycleIndex` do read -u6 { $analysisScript $inputPath $index $outputPath echo &>6 } done wait exec 6&>$-

執行

./job.sh ./analysis.py ./snap/snapshot_%03d.hdf5 3000 ./png/snap_

基本上執行完你就將數據文件轉化成圖片了.

接下來將圖片合併成視頻.

安裝ffmpeg,然後將png轉化成MP4

suod apt install ffmpeg cd cd Gadget-2.0.7/galaxy/png ffmpeg -i snap*.png ../movie.mp4

搞定.

說用Python的都死心了吧，速度絕對跟不上的。

乖乖的用C/Fortran還差不多。當年自己PC寫project跑拉格朗日點的龐加萊截面都要半個小時（大概是因為自己比較笨）。隔壁大氣我知道都是server跑Fortran來著。

用Python可以，前提是你要用既有的或者自己寫的C庫

去查一下n-body simulation

參考 EVE