N-Body simulation using the Particle Simulation Description Language

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 :)

--