monte carlo 1d radiative transfer code written in haskell using random-mersenne package and parallel computation technique

mc1drad_mersenne_parallel.hs

import System.Random.Mersenne

import Control.Parallel.Strategies hiding (parMap)

nt  = 4 -- number of thread
np  = 1000000 -- number of photon
tau = 1.0 -- cloud optical thickness
sal = 0.5 -- surface albedo
omg = 1.0 -- single scattering albedo
zu  = tau -- atmosphere upper boundary
zb  = 0.0 -- atmosphere lower boundary

main :: IO ()
main = do
    seed0 <- getStdGen
    resEval <- runEval $ parMap (world seed0) $ replicate nt $ fromIntegral $ np `div` nt
    (sumr, sumt) <- return $ mapTuple sum $ unzip $ resEval
    print $ [sumr / (fromIntegral np), sumt / (fromIntegral np)]

world :: MTGen -> Int -> IO (Double, Double)
world seed0 np' = loop 0 0.0 0.0 seed0
    where
        loop :: Int -> Double -> Double -> MTGen -> IO (Double, Double)
        loop t sumr sumt seed
            | t == np' = return (sumr, sumt)
            | otherwise = do
                (wr, wt) <- solver 1.0 zu (-1.0) 0.0 0.0 seed
                loop (t+1) (sumr + wr) (sumt + wt) seed

solver :: Double -> Double -> Double
          -> Double -> Double
          -> MTGen -> IO (Double, Double)
solver w z dir sumr sumt seed = do
    n <- random seed :: IO Double
    path <- return $ -log(n)
    nextphase (z + dir * path)
    where
        nextphase zn
            | zn <= zb = solver (w * sal) zb 1.0 sumr (sumt + w) seed
            | zn > zu = return (sumr + w, sumt)
            | otherwise  = do
                n2 <- random seed :: IO Double
                dnew <- return $ if n2 <= 0.5
                                   then 1.0 * dir
                                   else -1.0 * dir
                solver (w * omg) zn dnew sumr sumt seed

mapTuple :: (a -> b) -> (a, a) -> (b, b)
mapTuple f (a1, a2) = (f a1, f a2)

parMap :: (a -> IO b) -> [a] -> Eval (IO [b])
parMap _ [] = return $ return []
parMap f (a:as) = do
    b <- rpar $ f a
    bs <- parMap f as
    return $ do
        b' <- b
        bs' <- bs
        return (b':bs')

benchmark

for 10000000 photons

[rinx][~/mc1drad_mwc]% time ./mc1drad_mersenne_parallel +RTS -N2 -RTS
[0.6000630635562897,0.7998738728874206]
./mc1drad_mersenne_parallel +RTS -N2 -RTS  48.94s user 56.13s system 98% cpu 1:46.85 total
[rinx][~/mc1drad_mwc]% vi mc1drad_mersenne_parallel.hs
[rinx][~/mc1drad_mwc]% fg
fg: no current job
[rinx][~/mc1drad_mwc]% vi mc1drad_mersenne_parallel.hs
[rinx][~/mc1drad_mwc]% ghc -O3 -threaded mc1drad_mersenne_parallel.hs -o mc1drad_mersenne_parallel
[1 of 1] Compiling Main             ( mc1drad_mersenne_parallel.hs, mc1drad_mersenne_parallel.o )
Linking mc1drad_mersenne_parallel ...
[rinx][~/mc1drad_mwc]% time ./mc1drad_mersenne_parallel +RTS -N -RTS
[0.5999425475341796,0.8001149049316406]
./mc1drad_mersenne_parallel +RTS -N -RTS  53.07s user 87.47s system 153% cpu 1:31.34 total
[rinx][~/mc1drad_mwc]% fg
fg: no current job
[rinx][~/mc1drad_mwc]% vi mc1drad_mersenne_parallel.hs
[rinx][~/mc1drad_mwc]% ghc -O3 -threaded mc1drad_mersenne_parallel.hs -o mc1drad_mersenne_parallel
[1 of 1] Compiling Main             ( mc1drad_mersenne_parallel.hs, mc1drad_mersenne_parallel.o )
Linking mc1drad_mersenne_parallel ...
[rinx][~/mc1drad_mwc]% fg
fg: no current job
[rinx][~/mc1drad_mwc]% vi mc1drad_mersenne_parallel.hs
[rinx][~/mc1drad_mwc]% ghc -O3 -threaded mc1drad_mersenne_parallel.hs -o mc1drad_mersenne_parallel
[1 of 1] Compiling Main             ( mc1drad_mersenne_parallel.hs, mc1drad_mersenne_parallel.o )
Linking mc1drad_mersenne_parallel ...
[rinx][~/mc1drad_mwc]% time ./mc1drad_mersenne_parallel +RTS -N -RTS
[0.6000811416107178,0.7998377167785644]
./mc1drad_mersenne_parallel +RTS -N -RTS  48.24s user 86.78s system 153% cpu 1:27.91 total