On
the former entry (in Japanese), a free-fall simulation on point mass is shown as a starting point towards designing a language to describe particle simulation, or SPH.
On this entry, I introduce a improved Language and demonstrate N-Body Simulation on point mass using it.
Describing N-Body Simulation
N-Body simulation is represented in mathematical formula as below:
where xi, vi, ai and fi are position, velocity, acceleration and force on point mass i, respectively.
With the introduced language, as embedded DSL in Haskell, this simulation can be described as below:
nBody =
do define "x" (integral (var "v" <*> dt))
define "v" (integral (var "a" <*> dt))
defineWithType "a" (var "f" </> m) ExpTypeVector
define "f"
(let r = norm (var' "x" <-> var "x")
k = m <*> m <*> g </> r </> r
n = (var' "x" <-> var "x") </> r
in sigma (k <*> n))
initialConditionV "x" x0
where dt = constantS 0.01
m = constantS 1
g = constantS 9.8
x0 = V.fromList [ (0,0,0), (2,0,0), (0,2,0) ]
Explanation
define declaration defines an equation on physical amount on point mass i.
defineWithType declaration also defines an equation on physical amount, instructing type with it. When we use
define declaration, the type of the equation we define is inferred with other equations. On the other hand, if the type can not be inferred, we must instruct it using
defineWithType declaration. The types we can use are Scalar and Vector for now.
initialConditionV declaration is to give the initial condition of the physical amount defined with a "integral" expression. The last letter "V" means that the physical amount has Vector type. As the argument of it, an array of Vector is given as the initial condition of the position "x". Besides, the initial condition of the velocity "v" can be omitted, and an array of zero vector is given implicitly.
integral expression means integration, calculated explicitly inside the module.
var expression and
var' expression are variables indicated by the symbol given as its argument.
var' expression must be in a sigma expression. The difference in
var and
var' is described in the next paragraph.
sigma expression means summation of the calculated result of its operand expression on point mass j excluding itself. In the operand,
var "x" means xi and
var' "x" means xj.
dt,
m and
g are constant values lifted with the
constantS function. The last letter "S" shows that the constant has Scalar type.
Execution
The nBody defined in the code can be evaluated as below:
main = mapM_ (printRow . flip machineRegisterValueV "x")
$ take 100 $ runSimMachine nBody
* for the whole code, see
Simulation.hs in GitHub
About performance efficiency
For now, this language is implemented as interpreter, so the performance of codes written in it is not very fast. I've compared it with the similar N-Body simulation code written in plain Haskell (see
NBody.hs in GitHub). The result is below. Simulation is in the language and NBody in plain Haskell.
$ time ./Simulation > /dev/null
real 0m2.900s
user 0m2.869s
sys 0m0.020s
$ time ./NBody > /dev/null
real 0m0.846s
user 0m0.829s
sys 0m0.011s
The program in the language is x1/3 - x1/4 as fast as in plain Haskell.
To improve the performance, using Template Haskell, we can compile, in the broad sense, the language to plain Haskell and make it work more efficiently. Also, compiling it to generate C or CUDA code would be pretty nice :)
--
粒子法(SPH)のプログラムを解説したシリーズです。ソースコードも公開しています。
シミュレーションの結果をレンダリングして作った動画です。流体シミュレーションや剛体シミュレーションの動画を見ることができます。
