# How to make following CFD program more idiomatic?

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 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 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?

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.

• Thanks, I've misunderstood how where works in general. Also, I've ran hlint. Mar 15, 2012 at 14:51