3-body Problem in Haskell
So last week I implemented the n-body simulation in Erlang with message-passing, as an exploration of one-process-per-particle concurrency. Now I've got it up and running in Haskell, one version using STM and another using just MVars (individual atomic variables).
I found the shared-memory approach to be a more natural style than message-passing for this kind of physical model. My technique is simple: for each time step in simulation time, there is a shared vector of the states of the particles. Each process constantly tries to read the current moment values in sequence (blocking on those particles that are not filled in) and when it has read the whole thing, it performs its next-state computation on those states; it then writes its own next-state into the t+1 shared state vector. I want to reiterate that this technique works only because the processes necessarily work synchronously through simulation time.
Why MVars instead of STM? MVars implement a medium-power concurrency primitive—something like test-and-set—and it was plenty powerful enough to write this application without much fuss. Transactional memory is a significantly higher-level primitive, and its added flexibility wasn't necessary in this case. I'd like to make a loose argument that MVars are generally suffficient for the concurrency that's likely to be found in games.
The STM version works very similarly, but instead of just blocking on an unfilled particle-state, it retries. I believe that this is inefficient, because it's not necessary to repeat the reads that have already been performed at a given timestep: the variables in question are all single-assignment. Also, I suspect that the semantics of STM will lead to lots of retries, the threads tending to swarm around an unfilled variable and all of them retrying the whole transaction together. By contrast, the MVar implementation is efficient in the sense that when a state-variable is filled in, one waiting thread is woken, whose take/put action causes another waiting thread to be woken, etc.
Here's the code, for reference.
-- Ezra's hackneyed gravity sim inConcurrent Haskell, MVar version
import Control.Concurrent
(+:+) :: Num a => (a, a) -> (a, a) -> (a, a)
(+:+) (x1, y1) (x2, y2) = (x1 + x2, y1 + y2)
(*:*) :: Num a => a -> (a, a) -> (a, a)
(*:*) s (x, y) = (x * s, y * s)
force :: (Ord s, Fractional s) => (s, s) -> (s, s) -> (s, s)
force (ax, ay) (bx, by) = if distsq > 0 then (dx / distsq, dy / distsq) else (0, 0)
where dx = (bx - ax)
dy = (by - ay)
distsq = dx * dx + dy * dy
timestep :: (Ord s, Fractional s) => (s, s) -> (s, s) -> [(s, s)] -> ((s, s), (s, s))
timestep pos vel pos_vector = (pos +:+ vel, vel +:+ dv)
where dv = 0.001 *:* (foldr1 (+:+) (map (\pos' -> force pos pos') pos_vector))
readAllMVars :: [MVar a] -> IO [a]
readAllMVars vars = sequence (map readMVar vars)
rounds = 100
n = 3
new_pos_vec = sequence (map (\_-> newEmptyMVar) [1..n])
particle :: Int -> Int -> ((Float, Float), (Float, Float))
-> [MVar (Float, Float)] -> [[MVar (Float, Float)]]
-- -> ((Float, Float), (Float, Float))
-> IO ()
particle i t state currentPosns [] = return ()
particle i t (pos, vel) currentPosns (nextPosns:futurePosns) =
do {
nextState <- do {
posnsVec <- readAllMVars currentPosns;
(pos, vel) <- return(timestep pos vel posnsVec) ;
putMVar (nextPosns !! i) pos;
return (pos, vel)
};
particle i (t+1) nextState nextPosns futurePosns
}
-- end when futurePosns is exhausted?
glue g [] = ""
glue g [x] = x
glue g (h:t) = h ++ g ++ glue g t
watch [] = return ()
watch (state : futureStates) =
do { posns <- readAllMVars state;
putStr $ (glue ", " $ map show posns) ++ "\n";
watch futureStates
}
main =
-- the initial position and velocity of each particle. They all start
-- out at rest; their x co-ordinates are 1, 2, 3, 4, ... and their
-- y co-ords alternate +1, -1, +1, -1, ... This way they don't all
-- sit on a line, which would be boring.
let initialStates = [((fromIntegral proc, fromIntegral (proc `mod` 2) * 2 - 1), (0.00, 0.00)) | proc <- [0..n-1]] in
do {
-- stateVectors and initialPubStates are the shared variables where
-- the processes will communicate their state. initialPubStates is
-- the the initialized first position, whereas the stateVectors all
-- start out empty.
stateVectors <- sequence [new_pos_vec | t <- [1..rounds]] ;
initialPubStates <- sequence $ map (\(pos, vel) -> newMVar pos) initialStates;
sequence $ zipWith (\s i -> forkIO $ particle i 0 s initialPubStates stateVectors) initialStates [0..n-1];
watch stateVectors;
}
Comments
I must thank you for the efforts you've put in penning this website. I am hoping to check out the same high-grade blog posts by you later on as well. In truth, your creative writing abilities has inspired me to get my own, personal site now ;)
Posted by: john thomas financial | October 18, 2012 10:23 PM