Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

begin raising Haskell repa program from the dead #6

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
works but SLOW lol
raehik committed Feb 14, 2022
commit 83cbc1cd18292753f1b7f0f9cb6fbcdc994c1fe6
76 changes: 76 additions & 0 deletions haskell/repa-2022/app/Main.hs
Original file line number Diff line number Diff line change
@@ -0,0 +1,76 @@
{-# LANGUAGE NoMonomorphismRestriction #-}

module Main where

import Navier.Boundary
import Navier.Datadefs
import Navier.Init
import Navier.Params
import Navier.Output
import Navier.Simulation

import Data.Array.Repa

import System.Environment
import Debug.Trace
import Text.Printf

main = do
argv <- getArgs
let output_name = if ((length argv) >= 1) then argv!!0 else "foo"
output_freq = if ((length argv) >= 2) then read $ argv!!1 else 1
dims = Z :. (jmax+2) :. (imax+2)

-- initialise arrays
u = fromListUnboxed dims (Prelude.replicate ((imax+2)*(jmax+2)) ui)
v = fromListUnboxed dims (Prelude.replicate ((imax+2)*(jmax+2)) vi)
f = fromListUnboxed dims (Prelude.replicate ((imax+2)*(jmax+2)) 0.0)
g = fromListUnboxed dims (Prelude.replicate ((imax+2)*(jmax+2)) 0.0)
rhs = fromListUnboxed dims (Prelude.replicate ((imax+2)*(jmax+2)) 0.0)
p = fromListUnboxed dims (Prelude.replicate ((imax+2)*(jmax+2)) 0.0)

-- initialise the "flag" array, marking a circular peg
(flag, ibound) = init_flag imax jmax delx dely
ifluid = (imax * jmax) - ibound

-- initial boundary conditions on the velocity components
(u', v') = apply_boundary_conditions u v flag imax jmax ui vi
output = mainLoop 0.0 0 del_t 0.0 u' v' p f g rhs flag ifluid ibound output_freq

-- write frames out
write_ppm u v p flag imax jmax delx dely output_name 0 output_freq
mapM (\(u, v, p, i) ->
write_ppm u v p flag imax jmax delx dely output_name (i+1) output_freq) output

mainLoop
:: Double -> Int -> Double -> Double
-> Array U DIM2 Double -> Array U DIM2 Double
-> Array U DIM2 Double -> Array U DIM2 Double
-> Array U DIM2 Double -> Array U DIM2 Double
-> Array U DIM2 Int
-> Int -> Int -> Int
-> [(Array U DIM2 Double, Array U DIM2 Double, Array U DIM2 Double, Int)]
mainLoop t iters del_t res u v p f g rhs flag ifluid ibound output_freq
| t < t_end =
let
del_t' = set_timestamp_interval del_t u v
(f', g') = compute_tentative_velocity u v f g flag del_t'
rhs' = compute_rhs f' g' rhs flag del_t'

-- the heavy computational work mostly happens in poisson
(p', res', itersor) = poisson p rhs' flag ifluid res

(u', v') = update_velocity u v f' g' p' flag del_t'
(u'', v'') = apply_boundary_conditions u' v' flag imax jmax ui vi

-- output message
msg = printf "%d t:%.8f del_t:%.8f, SOR iters:%3d, res:%.8f, bcells:%d"
iters (t+del_t') del_t' itersor res' ibound
frameOutput = (u'', v'', p', iters)
in
trace msg $
if (iters `mod` output_freq == 0) then
frameOutput:(mainLoop (t+del_t') (iters+1) del_t' res' u'' v'' p' f' g' rhs' flag ifluid ibound output_freq)
else
mainLoop (t+del_t') (iters+1) del_t' res' u'' v'' p' f' g' rhs' flag ifluid ibound output_freq
| otherwise = [(u, v, p, iters)]
32 changes: 22 additions & 10 deletions haskell/repa-2022/package.yaml
Original file line number Diff line number Diff line change
@@ -1,16 +1,6 @@
name: repa2022
version: "0.1.0"

dependencies:
- base
- repa
- vector
- template-haskell
- directory

library:
source-dirs: src

# mostly Alexis King's 2018 recommended defaults
# (most can be replaced with GHC 9.2's GHC2021 language extension
default-extensions:
@@ -50,3 +40,25 @@ default-extensions:
- UndecidableInstances # honestly fine but...
- MagicHash # pretty much syntactic, but too weird
- ScopedTypeVariables # probs dangerous to have as default

dependencies:
- base
- repa
- vector
- template-haskell
- directory

library:
source-dirs: src

executables:
repa2022:
source-dirs: app
main: Main.hs
ghc-options:
- -Wall
- -threaded
- -rtsopts
- -with-rtsopts=-N
dependencies:
- repa2022
49 changes: 48 additions & 1 deletion haskell/repa-2022/repa2022.cabal
Original file line number Diff line number Diff line change
@@ -11,7 +11,6 @@ build-type: Simple
library
exposed-modules:
Common
Navier
Navier.Boundary
Navier.Datadefs
Navier.Init
@@ -61,3 +60,51 @@ library
, template-haskell
, vector
default-language: Haskell2010

executable repa2022
main-is: Main.hs
other-modules:
Paths_repa2022
hs-source-dirs:
app
default-extensions:
EmptyCase
FlexibleContexts
FlexibleInstances
InstanceSigs
MultiParamTypeClasses
PolyKinds
LambdaCase
DerivingStrategies
StandaloneDeriving
DeriveAnyClass
DeriveGeneric
DeriveDataTypeable
DeriveFunctor
DeriveFoldable
DeriveTraversable
DeriveLift
ImportQualifiedPost
StandaloneKindSignatures
DerivingVia
RoleAnnotations
TypeApplications
DataKinds
TypeFamilies
TypeOperators
BangPatterns
GADTs
DefaultSignatures
RankNTypes
UndecidableInstances
MagicHash
ScopedTypeVariables
ghc-options: -Wall -threaded -rtsopts -with-rtsopts=-N
build-depends:
base
, directory
, repa
, repa2022
, template-haskell
, vector
default-language: Haskell2010
13 changes: 2 additions & 11 deletions haskell/repa-2022/src/Common.hs
Original file line number Diff line number Diff line change
@@ -1,16 +1,7 @@
module Common where

import Data.Array.Repa
import Data.Vector.Unboxed

-- TODO likely a strict pair so.
data a :*: b = !a :*: !b
infixr 9 :*:

fsta :: a :*: b -> a
fsta (a :*: _) = a

snda :: a :*: b -> b
snda (_ :*: b) = b

toScalar :: Array U DIM0 Int -> Int
toScalar :: Source r a => Array r DIM0 a -> a
toScalar arr = index arr Z
72 changes: 0 additions & 72 deletions haskell/repa-2022/src/Navier.lhs

This file was deleted.

2 changes: 1 addition & 1 deletion haskell/repa-2022/src/Navier/Datadefs.hs
Original file line number Diff line number Diff line change
@@ -78,7 +78,7 @@ mkPat (l, x, y) = varP $ mkName l
mkPats vars (l, _, _) = Prelude.map (varP . mkName) (tie vars (repeat l))

inter [v] l = varP $ mkName (v++l)
inter (v:vs) l = infixP (varP $ mkName (v++l)) (mkName ":*:") (inter vs l)
inter (v:vs) l = infixP (varP $ mkName (v++l)) (mkName ",") (inter vs l)

--Prelude.map (varP . mkName) (tie vars (repeat l))

4 changes: 1 addition & 3 deletions haskell/repa-2022/src/Navier/Init.hs
Original file line number Diff line number Diff line change
@@ -4,6 +4,7 @@ module Navier.Init where

import Prelude hiding ( traverse )

import Common
import Data.Bits
import Data.Array.Repa
import Data.Array.Repa.Eval ( suspendedComputeP )
@@ -70,6 +71,3 @@ init_flag imax jmax delx dely =
then 1+y else y
total y x = y + x
in (flag''', ibound)
where
toScalar :: Array U DIM0 Int -> Int
toScalar arr = index arr Z
10 changes: 5 additions & 5 deletions haskell/repa-2022/src/Navier/Output.hs
Original file line number Diff line number Diff line change
@@ -21,12 +21,12 @@ calc_zeta
-> Int -> Int -> Double -> Double
-> ArrayU Double
calc_zeta uA vA flagA imax jmax delx dely =
computeS $ traverse (Data.Array.Repa.zipWith (:*:) flagA
(Data.Array.Repa.zipWith (:*:) uA vA)) id calcZeta
computeS $ traverse (Data.Array.Repa.zipWith (,) flagA
(Data.Array.Repa.zipWith (,) uA vA)) id calcZeta
where
flag = fsta
u = fsta . snda
v = snda . snda
flag = fst
u = fst . snd
v = snd . snd
calcZeta get c@(sh :. j :. i) =
if (inBounds i j && (i<=(imax-1)) && (j<=(jmax-1))) then
if (((flag $ get c) .&. _cf /= 0) &&
4 changes: 2 additions & 2 deletions haskell/repa-2022/src/Navier/Params.hs
Original file line number Diff line number Diff line change
@@ -36,9 +36,9 @@ inBounds i j = (i>0) && (i<(imax+1)) && (j>0) && (j<(jmax+1))
intCast :: Int -> Double
intCast = fromInteger . toInteger

ignoreBoundary :: Unbox a => Array U DIM2 a -> Array U DIM2 a
ignoreBoundary :: (Unbox a, Source r a) => Array r DIM2 a -> Array D DIM2 a
ignoreBoundary arr = let Z :. height :. width = extent arr
in computeS $ backpermute (Z :. (height-2) :. (width-2))
in backpermute (Z :. (height-2) :. (width-2))
(\(sh :. y :. x) -> (sh :. (y+1) :. (x+1))) arr

{-# INLINE obs #-}
Loading