2
\$\begingroup\$

Basically, I'm comparably new to Haskell, which means lacking both knowledge of advanced principles, available libraries and, maybe, some syntax.

For educational purposes I hacked together a code implementing Stam-style fluid dynamics solver. For now the code is running in constant space and linear time, which is exactly what's needed. But feeling remains that code I wrote is rather ugly.

I have serious concerns about whether this code is possible to read at all, but hope that somebody could point out most obvious idiosyncrasies.

import Data.Array.Unboxed
import Data.List
import Text.Printf

n = 20
timeStep = 0.1
numSteps = 100

indicies = [[0..n + 1], [1..n], [2..n-1]]

type Matrix = UArray (Int, Int) Float
data Bound = XHard | YHard | Soft
data State = State !Matrix !Matrix !Matrix

applyNTimes :: Int -> (a -> a) -> a -> a
applyNTimes times f val = foldl' (\m i -> f m) val [1..times]

showMatrix :: Matrix -> String
showMatrix m = foldl (++) "" [showRow m i ++ "\n" | i <- head indicies] where showRow m i = foldl (++) "" [printf "%.3f " (m ! (i, j)) | j <- head indicies]

showState :: State -> String
showState (State u v d) = showMatrix d

createMatrix :: [Int] -> ((Int, Int) -> Float) -> Matrix
createMatrix ind f = array ((0, 0), (n + 1, n + 1)) [((i, j), f (i, j)) | j <- ind, i <- ind] // [(ix, 0) | ix <- other] 
            where other = [(i, j) | i <- head indicies, j <- head indicies] \\ [(i, j) | i <- ind, j <- ind]

zeroMatrix :: Matrix
zeroMatrix = createMatrix (head indicies) (\ (i, j) -> 0.0)

add :: Matrix -> Matrix -> Matrix
add a b = createMatrix (head indicies) (add_cell a b) where add_cell a b ix = a ! ix + b ! ix

addPecularSource :: State -> State
addPecularSource (State u v d) = State u v (d `add` createMatrix (head indicies) spike) where spike ix = if ix == (n `div` 2, n `div` 2) then 0.1 else 0.0

advance :: State -> State
advance s = solveDensity . addPecularSource . solveVelocity $ s

main :: IO()
main = putStrLn . showState . applyNTimes numSteps advance $ State z z z where z = zeroMatrix

solveDensity :: State -> State
solveDensity (State u v d) = State u v d' where d' = advect Soft (diffuse d diffusion) u v where diffusion = 0.0 :: Float

setCorners :: Matrix -> Matrix
setCorners a = a // [((n + 1, n + 1),   (a ! (n , n + 1) + a ! (n + 1, n)) / 2), 
             ((n + 1, 0),   (a ! (n , 0) + a ! (n + 1, 1)) / 2),
             ((0, n + 1),   (a ! (0, n) + a ! (1, n + 1)) / 2),
             ((0, 0),   (a ! (1 , 0) + a ! (0, 1)) / 2)]

diffuse :: Matrix -> Float -> Matrix
diffuse c0 diff = linearSolver c0 a (1 + 4 * a) where a = timeStep * diff * fromIntegral n * fromIntegral n

setBoundary :: Bound -> Matrix -> Matrix
setBoundary b m = setCorners (setSides b m)

setSides :: Bound -> Matrix -> Matrix
setSides XHard m = setSides Soft m // ([((0, i), - m ! (0, i)) | i <- indicies !! 1] ++ [((n + 1, i), - m ! (n + 1, i)) | i <- indicies !! 1])
setSides YHard m = setSides Soft m // ([((i, 0), - m ! (i, 0)) | i <- indicies !! 1] ++ [((i, n + 1), - m ! (i, n + 1)) | i <- indicies !! 1])
setSides Soft m = m // ([((0, i), m ! (1, i)) | i <- indicies !! 1] ++ [((n + 1, i), m ! (n, i)) | i <- indicies !! 1] ++
            [((i, 0), m ! (i, 1)) | i <- indicies !! 1] ++ [((i, n + 1), m ! (i, n)) | i <- indicies !! 1])

linearSolver :: Matrix -> Float -> Float -> Matrix
linearSolver x0 a c = applyNTimes 20 (iterFun x0 a c) zeroMatrix 
            where iterFun x0 a c m = setBoundary Soft (createMatrix (indicies !! 1) (update m x0 a c) ) 
                where update m x0 a c (i, j) = (a * (m ! (i - 1, j) + m ! (i + 1, j) + m ! (i, j - 1) + m ! (i, j + 1)) + x0 ! (i, j)) / c

curl :: Matrix -> Matrix  -> Matrix
curl u v = createMatrix (indicies !! 1) gen where gen (i, j) = (u ! (i, j + 1) - u ! (i, j - 1) - v ! (i + 1, j) + v ! (i - 1, j)) / 2

archCell :: Matrix -> Float -> Float -> (Int, Int) -> Float
archCell d gravity push (i, j) = gravity * d ! (i, j) + push * (d ! (i, j) - ambient d)
                    where ambient d = sum [d ! (i, j) | i <- indicies !! 1, j <- indicies !! 1]

archimedes :: Matrix -> Matrix
archimedes d = createMatrix (head indicies) (archCell d gravity push)
        where (gravity, push) = (0.001 :: Float, 0.005 :: Float)

solveVelocity :: State -> State
solveVelocity (State u v d) = State u' v' d
                where (u', v') = project (advect YHard uProjected uProjected vProjected) (advect XHard vProjected uProjected vProjected)
                    where (uProjected, vProjected) = project uDiffused vDiffused
                        where (uDiffused, vDiffused) = (diffuse uVorticityConfined viscosity, diffuse (vVorticityConfined `add` archimedes d) viscosity)
                            where ((uVorticityConfined, vVorticityConfined), viscosity) = (vorticityConfimentForce u v, 0.0)

advectCell :: Matrix -> Matrix -> Matrix -> (Int, Int) -> Float
advectCell d0 du dv (i, j) = s0 * t0 * d0 ! (floor x, floor y) + 
               s0 * t1 * d0 ! (floor x, floor y + 1) + 
               s1 * t0 * d0 ! (floor x + 1, floor y) + 
               s1 * t1 * d0 ! (floor x + 1, floor y + 1)
                where (s0, t0, s1, t1, x, y) = (1.0 - x' + xx, 1.0 - y' + yy, x' - xx, y' - yy, x', y')
                    where (x', y', xx, yy) = (x, y, fromIntegral (floor x), fromIntegral (floor y))
                        where (x, y) = (sort [0.5 + fromIntegral n, 0.5, fromIntegral i - dt0 * du ! (i, j)] !! 1,
                                        sort [0.5 + fromIntegral n, 0.5, fromIntegral j - dt0 * dv ! (i, j)] !! 1)
                                where dt0 = timeStep * fromIntegral n

advect :: Bound -> Matrix -> Matrix -> Matrix -> Matrix
advect b d0 du dv = setBoundary b (createMatrix (indicies !! 1) (advectCell d0 du dv))

projectCellX :: Matrix -> Matrix -> (Int, Int) -> Float
projectCellX x p (i, j) = x ! (i, j) - 0.5 * fromIntegral n * (p!(i + 1, j) - p!(i - 1, j))

projectCellY :: Matrix -> Matrix -> (Int, Int) -> Float
projectCellY y p (i, j) = y ! (i, j) - 0.5 * fromIntegral n * (p!(i, j + 1) - p!(i, j - 1))

project :: Matrix -> Matrix -> (Matrix, Matrix)
project x y = (setBoundary YHard x', setBoundary XHard y')
        where (x', y') = (createMatrix (indicies !! 1) (projectCellX x p), createMatrix (indicies !! 1) (projectCellY y p))
            where p = linearSolver (setBoundary Soft diverg) 1 4
                  where diverg = createMatrix (indicies !! 1) (computeDiv x y)
                             where computeDiv x y (i, j) = -0.5 * (x ! (i + 1, j) - x ! (i - 1, j) + y ! (i, j + 1) - y ! (i, j - 1)) / fromIntegral n

dwdx :: Matrix -> (Int, Int) -> Float
dwdx c (i, j) = (abs (c ! (i + 1, j)) - abs (c ! (i - 1, j))) / 2

dwdy :: Matrix -> (Int, Int) -> Float
dwdy c (i, j) = (abs (c ! (i, j + 1)) - abs (c ! (i, j - 1))) / 2

len :: Matrix -> (Int, Int) -> Float
len c (i, j) = sqrt(dwdx c (i, j) * dwdx c (i, j) + dwdy c (i, j) * dwdy c (i, j)) + eps where eps = 1e-6

funX :: Matrix -> (Int, Int) -> Float
funX c (i, j) = - c ! (i, j) * dwdy c (i, j) / len c (i, j)

funY :: Matrix -> (Int, Int) -> Float
funY c (i, j) = c ! (i, j) * dwdx c (i, j) / len c (i, j)

vorticityConfimentForce :: Matrix -> Matrix -> (Matrix, Matrix)
vorticityConfimentForce u v = (createMatrix (indicies !! 2) (funX (curl u v)), createMatrix (indicies !! 2) (funY (curl u v)))

GHC 6.12.3 was used. Bang-patterns were added via -XBangPatterns.

My question is: how to change this code to be more comprehensible and idiomatic?

\$\endgroup\$

migrated from stackoverflow.com Mar 14 '12 at 11:55

This question came from our site for professional and enthusiast programmers.

2
\$\begingroup\$

Get rid of all the cascading where clauses, e.g.

project x y = (setBoundary YHard x', setBoundary XHard y')
        where (x', y') = (createMatrix (indicies !! 1) (projectCellX x p), createMatrix (indicies !! 1) (projectCellY y p))
           where p = linearSolver (setBoundary Soft diverg) 1 4
                  where diverg = createMatrix (indicies !! 1) (computeDiv x y)
                             where computeDiv x y (i, j) = -0.5 * (x ! (i + 1, j) - x ! (i - 1, j) + y ! (i, j + 1) - y ! (i, j - 1)) / fromIntegral n

All of those can be in a single where block, and every declaration in a let or where block are in scope in the other declarations in that block.

Then I would run hlint on your code and take a look at what it suggests.

\$\endgroup\$
  • \$\begingroup\$ Thanks, I've misunderstood how where works in general. Also, I've ran hlint. \$\endgroup\$ – user1268256 Mar 15 '12 at 14:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.