From d9f02b03dcbe77e6324a5e54a54ba8879a8713f5 Mon Sep 17 00:00:00 2001 From: Ben Orchard Date: Thu, 10 Feb 2022 17:10:00 +0000 Subject: [PATCH 1/3] begin raising Haskell repa program from the dead Some concepts have been moved around the Haskell ecosystem, and Repa has had a bunch of breaking changes that need working through. Partway through the Repa changes. Might need need to rewrite it all to be monadic due to the changes... --- haskell/repa-2022/.gitignore | 1 + haskell/repa-2022/README.md | 9 + haskell/repa-2022/package.yaml | 52 +++ haskell/repa-2022/repa2022.cabal | 62 ++++ haskell/repa-2022/src/Navier.lhs | 72 ++++ haskell/repa-2022/src/Navier/Boundary.lhs | 122 +++++++ haskell/repa-2022/src/Navier/Datadefs.hs | 125 +++++++ haskell/repa-2022/src/Navier/Init.hs | 71 ++++ haskell/repa-2022/src/Navier/Output.hs | 62 ++++ haskell/repa-2022/src/Navier/Params.hs | 58 ++++ haskell/repa-2022/src/Navier/Simulation.lhs | 351 ++++++++++++++++++++ haskell/repa-2022/stack.yaml | 5 + haskell/repa-2022/stack.yaml.lock | 19 ++ 13 files changed, 1009 insertions(+) create mode 100644 haskell/repa-2022/.gitignore create mode 100644 haskell/repa-2022/README.md create mode 100644 haskell/repa-2022/package.yaml create mode 100644 haskell/repa-2022/repa2022.cabal create mode 100644 haskell/repa-2022/src/Navier.lhs create mode 100644 haskell/repa-2022/src/Navier/Boundary.lhs create mode 100644 haskell/repa-2022/src/Navier/Datadefs.hs create mode 100644 haskell/repa-2022/src/Navier/Init.hs create mode 100644 haskell/repa-2022/src/Navier/Output.hs create mode 100644 haskell/repa-2022/src/Navier/Params.hs create mode 100644 haskell/repa-2022/src/Navier/Simulation.lhs create mode 100644 haskell/repa-2022/stack.yaml create mode 100644 haskell/repa-2022/stack.yaml.lock diff --git a/haskell/repa-2022/.gitignore b/haskell/repa-2022/.gitignore new file mode 100644 index 0000000..6fabf46 --- /dev/null +++ b/haskell/repa-2022/.gitignore @@ -0,0 +1 @@ +/.stack-work/ diff --git a/haskell/repa-2022/README.md b/haskell/repa-2022/README.md new file mode 100644 index 0000000..c1230e4 --- /dev/null +++ b/haskell/repa-2022/README.md @@ -0,0 +1,9 @@ + * `computeP` places in a monad??? But any monad????? Don't get. Paper lies. + Using `computeS` instead, but less efficient. + * Think about where to force (compute), what reps combinators should use (e.g. + `Params.ignoreBoundary`) + * Similarly, had to use `foldS` because `foldP` requires a monad. Why doesn't + the upgrade paper describe the impacts of this change, or how to recover + original behaviour? + * I think I have to rewrite the *whole* thing to be parametric in some + monad... diff --git a/haskell/repa-2022/package.yaml b/haskell/repa-2022/package.yaml new file mode 100644 index 0000000..a9c8c94 --- /dev/null +++ b/haskell/repa-2022/package.yaml @@ -0,0 +1,52 @@ +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: +# essential +- EmptyCase +- FlexibleContexts +- FlexibleInstances +- InstanceSigs +- MultiParamTypeClasses +- PolyKinds +- LambdaCase +# deriving-related +- DerivingStrategies +- StandaloneDeriving +- DeriveAnyClass +- DeriveGeneric +- DeriveDataTypeable +- DeriveFunctor +- DeriveFoldable +- DeriveTraversable +- DeriveLift +# essential syntax but too recent +- ImportQualifiedPost # 8.10 +- StandaloneKindSignatures # 8.10 +- DerivingVia # 8.6 +# less essential but still gimmes +- RoleAnnotations +- TypeApplications +- DataKinds +- TypeFamilies +- TypeOperators +- BangPatterns +- GADTs +- DefaultSignatures +- RankNTypes +# extra +- UndecidableInstances # honestly fine but... +- MagicHash # pretty much syntactic, but too weird +- ScopedTypeVariables # probs dangerous to have as default diff --git a/haskell/repa-2022/repa2022.cabal b/haskell/repa-2022/repa2022.cabal new file mode 100644 index 0000000..82ddb99 --- /dev/null +++ b/haskell/repa-2022/repa2022.cabal @@ -0,0 +1,62 @@ +cabal-version: 1.12 + +-- This file has been generated from package.yaml by hpack version 0.34.4. +-- +-- see: https://github.com/sol/hpack + +name: repa2022 +version: 0.1.0 +build-type: Simple + +library + exposed-modules: + Navier + Navier.Boundary + Navier.Datadefs + Navier.Init + Navier.Output + Navier.Params + Navier.Simulation + other-modules: + Paths_repa2022 + hs-source-dirs: + src + 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 + build-depends: + base + , directory + , repa + , template-haskell + , vector + default-language: Haskell2010 diff --git a/haskell/repa-2022/src/Navier.lhs b/haskell/repa-2022/src/Navier.lhs new file mode 100644 index 0000000..03fcf85 --- /dev/null +++ b/haskell/repa-2022/src/Navier.lhs @@ -0,0 +1,72 @@ +> {-# LANGUAGE NoMonomorphismRestriction #-} + +> module Navier 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 = fromList dims (Prelude.replicate ((imax+2)*(jmax+2)) ui) +> v = fromList dims (Prelude.replicate ((imax+2)*(jmax+2)) vi) +> f = fromList dims (Prelude.replicate ((imax+2)*(jmax+2)) 0.0) +> g = fromList dims (Prelude.replicate ((imax+2)*(jmax+2)) 0.0) +> rhs = fromList dims (Prelude.replicate ((imax+2)*(jmax+2)) 0.0) +> p = fromList 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 DIM2 Double -> Array DIM2 Double +> -> Array DIM2 Double -> Array DIM2 Double -> Array DIM2 Double -> Array DIM2 Double +> -> Array DIM2 Int -> Int -> Int -> Int +> -> [(Array DIM2 Double, Array DIM2 Double, Array 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)] diff --git a/haskell/repa-2022/src/Navier/Boundary.lhs b/haskell/repa-2022/src/Navier/Boundary.lhs new file mode 100644 index 0000000..f838d85 --- /dev/null +++ b/haskell/repa-2022/src/Navier/Boundary.lhs @@ -0,0 +1,122 @@ +> {-# LANGUAGE NoMonomorphismRestriction #-} +> {-# LANGUAGE QuasiQuotes #-} + +> module Navier.Boundary where + +> import Data.Bits +> import Data.Array.Repa + +> import Navier.Datadefs +> import Navier.Params + +> import Text.Printf +> import Debug.Trace + +> apply_boundary_conditions :: Array DIM2 Double -> Array DIM2 Double -> +> Array DIM2 Int -> Int -> Int -> +> Double -> Double -> +> (Array DIM2 Double, Array DIM2 Double) +> apply_boundary_conditions (u@Manifest{}) (v@Manifest{}) +> (flag@Manifest{}) imax jmax ui vi = +> let ui'@Manifest{} = force $ traverse u id update +> where update get c@(sh :. j :. i) +> -- Fluid freely flows in from the west +> | (i==0) = get (sh :. j :. 1) +> -- Fluid freely flows out to the east +> | (i==imax) = get (sh :. j :. (imax-1)) +> | otherwise = get c +> u'@Manifest{} = force $ traverse ui' id update +> where update get c@(sh :. j :. i) +> -- fluid flows freely in the horizontal direction +> | (j==(jmax+1)) = get (sh :. jmax :. i) +> | (j==0) = get (sh :. 1 :. i) +> | otherwise = get c +> vi'@Manifest{} = force $ traverse v id update +> where update get c@(sh :. j :. i) +> -- Fluid freely flows in from the west +> | (i==0) = get (sh :. j :. 1) +> -- Fluid freely flows out to the east +> | (i==(imax+1)) = get (sh :. j :. imax) +> | otherwise = get c +> v'@Manifest{} = force $ traverse vi' id update +> where update get c@(sh :. j :. i) +> -- vertical velocity approaches 0 at the north and south +> | (j==jmax) = 0.0 +> | (j==0) = 0.0 +> | otherwise = get c +> -- Apply no-slip boundary conditions to cells that are adjacent to +> -- internal obstacle cells. This forces the u and v velocity to +> -- tend towards zero in these cells. +> u''@Manifest{} = force $ traverse2 u' flag (\x y -> x) update +> where +> update getU getFlag c@(sh :. j :. i) = +> if (inBounds i j && (getFlag c .&. _bnsew /= 0)) then +> case (getFlag c) of +> [p|_bn|] -> -getU (sh :. (j+1) :. i) +> [p|_be|] -> 0.0 +> [p|_bne|] -> 0.0 +> [p|_bse|] -> 0.0 +> [p|_bnw|] -> -getU (sh :. (j+1) :. i) +> [p|_bs|] -> -getU (sh :. (j-1) :. i) +> [p|_bsw|] -> -getU (sh :. (j-1) :. i) +> _ -> getU c +> else +> getU c +> u'''@Manifest{} = force $ traverse2 u'' flag (\x y -> x) update +> where +> update getU getFlag c@(sh :. j :. i) = +> if (((inBounds i j || i==0) && (i<=(imax-1))) && +> (getFlag (sh :. j :. (i+1)) .&. _bnsew /= 0)) then +> case (getFlag (sh :. j :. (i+1))) of +> [p|_bn|] -> -getU (sh :. (j+1) :. i) +> [p|_bw|] -> 0.0 +> [p|_bne|] -> -getU (sh :. (j+1) :. i) +> [p|_bsw|] -> 0.0 +> [p|_bnw|] -> 0.0 +> [p|_bs|] -> -getU (sh :. (j-1) :. i) +> [p|_bse|] -> -getU (sh :. (j-1) :. i) +> _ -> getU c +> else +> getU c +> v''@Manifest{} = force $ traverse2 v' flag (\x y -> x) update +> where +> update getV getFlag c@(sh :. j :. i) = +> if (inBounds i j && (getFlag c .&. _bnsew /= 0)) then +> case (getFlag c) of +> [p|_bn|] -> 0.0 +> [p|_be|] -> -getV (sh :. j :. (i+1)) +> [p|_bne|] -> 0.0 +> [p|_bse|] -> -getV (sh :. j :. (i+1)) +> [p|_bnw|] -> 0.0 +> [p|_bw|] -> -getV (sh :. j :. (i-1)) +> [p|_bsw|] -> -getV (sh :. j :. (i-1)) +> _ -> getV c +> else +> getV c +> v'''@Manifest{} = force $ traverse2 v'' flag (\x y -> x) update +> where +> update getV getFlag c@(sh :. j :. i) = +> if (((inBounds i j || j==0) && j<=(jmax-1)) && +> (getFlag (sh :. (j+1) :. i) .&. _bnsew /= 0)) then +> case (getFlag (sh :. (j+1) :. i)) of +> [p|_be|] -> -getV (sh :. j :. (i+1)) +> [p|_bs|] -> 0.0 +> [p|_bne|] -> -getV (sh :. j :. (i+1)) +> [p|_bse|] -> 0.0 +> [p|_bsw|] -> 0.0 +> [p|_bw|] -> -getV (sh :. j :. (i-1)) +> [p|_bnw|] -> -getV (sh :. j :. (i-1)) +> _ -> getV c +> else +> getV c +> -- Finally, fix the horizontal velocity at the western edge to have +> -- a continual flow of fluid into the simulation. +> u4@Manifest{} = force $ traverse u''' id update +> where update get c@(sh :. j :. i) +> | (i==0 && j>=1 && j<=jmax) = ui +> | otherwise = get c +> v4@Manifest{} = force $ traverse v''' id update +> where update get c@(sh :. j :. i) +> | (i==0 && j<=jmax) = 2*vi-(get (sh :. j :. 1)) +> | otherwise = get c +> in (u4, v4) \ No newline at end of file diff --git a/haskell/repa-2022/src/Navier/Datadefs.hs b/haskell/repa-2022/src/Navier/Datadefs.hs new file mode 100644 index 0000000..fd25596 --- /dev/null +++ b/haskell/repa-2022/src/Navier/Datadefs.hs @@ -0,0 +1,125 @@ +{-# LANGUAGE NoMonomorphismRestriction #-} +{-# LANGUAGE TemplateHaskell #-} + +module Navier.Datadefs where + +import Prelude hiding ( (++) ) + +import Data.Bits +import Language.Haskell.TH +import Language.Haskell.TH.Quote +import Data.Array.Repa + +import Debug.Trace + +-- 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 + +_cb = 0x0000 +_bn = 0x0001 +_bs = 0x0002 +_bw = 0x0004 +_be = 0x0008 +_bnw = _bn .|. _bw +_bsw = _bs .|. _bw +_bne = _bn .|. _be +_bse = _bs .|. _be +_bnsew = _bn .|. _bs .|. _be .|. _bw +_cf = 0x0010 + +p :: QuasiQuoter +p = QuasiQuoter undefined constPat undefined undefined + +constPat :: String -> PatQ +constPat varName = litP $ integerL $ case varName of + "_be" -> _be + "_bn" -> _bn + "_bs" -> _bs + "_bw" -> _bw + "_bne" -> _bne + "_bse" -> _bse + "_bsw" -> _bsw + "_bnw" -> _bnw + _ -> trace "FAIL!" 0 + +{- + +stencil :: QuasiQuoter +stencil = QuasiQuoter stencilExp stencilPat undefined stencilDecl + +tie _ [] = [] +tie [] _ = [] +tie (x:xs) (y:ys) = (x++y):(tie xs ys) + +pos :: [(String, Int, Int)] +pos = [("_tl", -1, -1), ("_t", 0, -1), ("_tr", 1, -1), + ("_l", -1, 0), ("_c", 0, 0), ("_r", 1, 0), + ("_bl", -1, 1), ("_b", 0, 1), ("_br", 1, 1)] + +stencilDecl :: String -> Q [Dec] +stencilDecl varsString = + do vars <- return $ words varsString + mapM (binder vars) pos + +binder vars (l, x, y) = + let pat = Prelude.map (varP . mkName) (tie vars (repeat l)) + exp = appE (varE $ mkName "get") + [| ($(varE $ mkName "sh")) :. + ($(varE $ mkName "j") + y) :. + ($(varE $ mkName "i") + x) |] + in valD (tupP pat) (normalB exp) [] + +stencilPat :: String -> Q Pat +stencilPat varsString + = do vars <- return $ words varsString + tupP [tupP (Prelude.map (mkPat) pos), + tupP (Prelude.map (varP . mkName) vars), + tupP (Prelude.concatMap (mkPats vars) pos)] + +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) + +--Prelude.map (varP . mkName) (tie vars (repeat l)) + +stencilExp :: String -> Q Exp +stencilExp varsString = do vars <- return $ words varsString + tupE [tupE (Prelude.map mkExps pos), + tupE (accessors vars), + tupE (Prelude.concatMap (mappings vars) pos)] + +mkExps (_, x, y) = + appE (varE $ mkName "get") + [| ($(varE $ mkName "sh")) :. + ($(varE $ mkName "j") + y) :. + ($(varE $ mkName "i") + x) |] + +accessors vars = Prelude.map conv (accs vars) + +--(iterate add1 [False])) + +mappings vars (l, x, y) = Prelude.map (\v -> appE (varE $ mkName v) + (varE $ mkName l)) vars + +conv xs = Prelude.foldr1 (\x y -> infixE (Just x) (varE $ mkName ".") (Just y)) + (Prelude.map (varE . mkName) xs) + +add1 :: [Bool] -> [Bool] +add1 [] = [True] +add1 (x:xs) = if x then (not x):(add1 xs) else True:xs + +accs vars = accs' (tail vars) ["fsta"] +accs' [x] acc = [acc, "snda":(tail acc)] +accs' (x:xs) acc = acc:(accs' xs (acc++["snda"])) + +-} diff --git a/haskell/repa-2022/src/Navier/Init.hs b/haskell/repa-2022/src/Navier/Init.hs new file mode 100644 index 0000000..d156b7c --- /dev/null +++ b/haskell/repa-2022/src/Navier/Init.hs @@ -0,0 +1,71 @@ +{-# LANGUAGE NoMonomorphismRestriction #-} + +module Navier.Init where + +import Prelude hiding ( traverse ) + +import Data.Bits +import Data.Array.Repa + +import Navier.Datadefs +import Navier.Params + +init_flag :: Int -> Int -> Double -> Double -> (Array U DIM2 Int, Int) +init_flag imax jmax delx dely = + let + dims = Z :. (jmax+2) :. (imax+2) + flag = suspendedComputeP $ fromList dims (Prelude.replicate ((imax+2)*(jmax+2)) (0::Int)) + -- Make a circular obstacle + mx = 20.0/41.0*(intCast jmax)*dely + my = mx + rad1 = 5.0/41.0*(intCast jmax)*dely + flag' = suspendedComputeP $ traverse flag id update + where + update get c@(sh :. j :. i) = + if (inBounds i j) then + let x = ((intCast i)-0.5)*delx - mx + y = ((intCast j)-0.5)*dely - my + in + if (x*x + y*y <= rad1*rad1) then _cb else _cf + else + get c + -- Make the north, sourth, east, and west boundary cells + flagi' = suspendedComputeP $ traverse flag' id update + where + update get c@(sh :. j :. i) = + if (j==0) then _cb + else if (j==(jmax+1)) then _cb + else get c + flag'' = suspendedComputeP $ traverse flagi' id update + where + update get c@(sh :. j :. i) = + if (j>=1 && j<=jmax) then + if (i==0) then _cb + else if (i==(imax+1)) then _cb + else get c + else + get c + -- Flag boundaries near obstacles + flag''' = suspendedComputeP $ traverse flag'' id update + where + update get c@(sh :. j :. i) = + if (inBounds i j && (get c .&. _cf) == 0) then + let + l = if (get (sh :. j :. (i-1)) .&. _cf /= 0) + then _bw else 0 + r = if (get (sh :. j :. (i+1)) .&. _cf /= 0) + then _be else 0 + b = if (get (sh :. (j-1) :. i) .&. _cf /= 0) + then _bs else 0 + t = if (get (sh :. (j+1) :. i) .&. _cf /= 0) + then _bn else 0 + in + (get c) .|. l .|. r .|. b .|. t + else + get (sh :. j :. i) + -- Count boundary cells + ibound = toScalar (foldS total 0 (foldS count 0 (ignoreBoundary flag'''))) + where count y x = if (x .&. _cf == 0) + then 1+y else y + total y x = y + x + in (flag''', ibound) diff --git a/haskell/repa-2022/src/Navier/Output.hs b/haskell/repa-2022/src/Navier/Output.hs new file mode 100644 index 0000000..69e466f --- /dev/null +++ b/haskell/repa-2022/src/Navier/Output.hs @@ -0,0 +1,62 @@ +module Navier.Output where + +import Prelude hiding ( traverse ) + +import Data.Bits +import Data.Char +import Data.Array.Repa + +import System.Directory +import System.IO +import Text.Printf + +import Navier.Datadefs +import Navier.Params + +type ArrayU a = Array U DIM2 a + +calc_zeta + :: ArrayU Double -> ArrayU Double -> ArrayU Int + -> 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 + where + flag = fsta + u = fsta . snda + v = snda . snda + calcZeta get c@(sh :. j :. i) = + if (inBounds i j && (i<=(imax-1)) && (j<=(jmax-1))) then + if (((flag $ get c) .&. _cf /= 0) && + ((flag $ get (sh :. j :. (i+1))) .&. _cf /= 0) && + ((flag $ get (sh :. (j+1) :. i)) .&. _cf /= 0) && + ((flag $ get (sh :. (j+1) :. (i+1))) .&. _cf /= 0)) then + ((u $ get (sh :. (j+1) :. i)) - (u $ get c))/dely + - ((v $ get (sh :. j :. (i+1))) - (v $ get c))/delx + else + 0.0 + else + 0.0 + +write_ppm + :: ArrayU Double -> ArrayU Double -> ArrayU Double -> ArrayU Int + -> Int -> Int -> Double -> Double -> String -> Int -> Int + -> IO () +write_ppm u v p flag imax jmax delx dely outname iters freq = do + createDirectoryIfMissing True outname + let filename = printf "%s/%06d.ppm" outname (iters `div` freq) + file <- openBinaryFile filename WriteMode + let zeta = calc_zeta u v flag imax jmax delx dely + hPutStr file $ "P6 " <> (show imax) <> " " <> (show jmax) <> " 255\n" + mapM (\(f, z) -> if (f .&. _cf == 0) then writeRGB file 0 255 0 + else writeBW file (((abs (z/12.6))**0.4)*255)) + (zip (toList $ ignoreBoundary flag) + (toList $ ignoreBoundary zeta)) + hClose file + +writeBW file b = writeRGB file b b b +writeRGB file r g b = do hPutChar file $ chr $ truncate r + hPutChar file $ chr $ truncate g + hPutChar file $ chr $ truncate b + diff --git a/haskell/repa-2022/src/Navier/Params.hs b/haskell/repa-2022/src/Navier/Params.hs new file mode 100644 index 0000000..2a3c6ee --- /dev/null +++ b/haskell/repa-2022/src/Navier/Params.hs @@ -0,0 +1,58 @@ +module Navier.Params where + +import Prelude hiding ( traverse ) + +import Data.Bits +import Data.Array.Repa +import Navier.Datadefs +import Data.Vector.Unboxed.Base ( Unbox ) + +-- boundary flags +xlength = 22.0 -- Width of simulated domain +ylength = 4.1 -- Height of simulated domain +imax = 660 -- Number of cells horizontally +jmax = 120 -- Number of cells vertically + +delx = xlength/660.0 +dely = ylength/120.0 + +t_end = 40.0 -- Simultion runtime +del_t = 0.003 -- (default) duration of each timestep +tau = 0.5 -- safety factor for timestep control + +itermax = 10::Int -- max number of SOR iterations +eps = 0.001 -- stopping error threshold for SOR +omega = 1.7 -- relaxation parameter for SOR +gamma = 0.9 -- upwind differencing factor in PDE discretisation + +re = 150.0 -- Reynolds number +ui = 1.0 -- Initial X velocity (West-east flow) +vi = 0.0 -- Initial Y velocity + +{-# INLINE inBounds #-} +inBounds :: Int -> Int -> Bool +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 arr = let Z :. height :. width = extent arr + in computeS $ backpermute (Z :. (height-2) :. (width-2)) + (\(sh :. y :. x) -> (sh :. (y+1) :. (x+1))) arr + +{-# INLINE obs #-} +obs x = if (x .&. _cf /= 0) then 1 else 0 + +{-# INLINE obstacle #-} +obstacle acc get i = if ((acc $ get i) .&. _cf /= 0) then 1 else 0 + +checksum_char :: (Source r Int, Monad m) => Array r DIM2 Int -> m Int +checksum_char arr = sumAllP $ traverse arr id update + where update get c@(sh :. j :. i) = + (get c)*(i) + +checksum :: (Source r Double, Monad m) => Array r DIM2 Double -> m Double +checksum arr = sumAllP $ traverse arr id update + where update get c@(sh :. j :. i) = + (get c)*(intCast $ (i*jmax)+j) diff --git a/haskell/repa-2022/src/Navier/Simulation.lhs b/haskell/repa-2022/src/Navier/Simulation.lhs new file mode 100644 index 0000000..5ccd4c3 --- /dev/null +++ b/haskell/repa-2022/src/Navier/Simulation.lhs @@ -0,0 +1,351 @@ + {-# LANGUAGE QuasiQuotes #-} + {-# LANGUAGE TypeOperators #-} + +> module Navier.Simulation where + +> import Navier.Datadefs +> import Navier.Params +> import Data.Array.Repa + +> import Debug.Trace +> import Data.Bits + +> import Text.Printf + +> import Data.Array.Parallel.Base ((:*:)(..)) + +> set_timestamp_interval del_t (u@Manifest{}) (v@Manifest{}) = +> -- del_t satisfying CFL conditions +> -- else no time stepsize control +> if (tau >= 1.0e-10) then +> let +> umax = foldAll (\x y -> max (abs x) (abs y)) (1.0e-10) +> (backpermute (Z :. (jmax-1) :. imax) +> (\(sh :. j :. i) -> (sh :. (j+1) :. i)) u) +> vmax = foldAll (\x y -> max (abs x) (abs y)) (1.0e-10) +> (backpermute (Z :. jmax :. (imax-1)) +> (\(sh :. j :. i) -> (sh :. j :. (i+1))) v) +> deltu = delx/umax +> deltv = dely/vmax +> deltre = 1.0/(1.0/(delx*delx)+1/(dely*dely))*re/2.0 +> in +> tau * (if (deltu else min deltv deltre) +> else +> del_t + +> compute_tentative_velocity (uA@Manifest{}) (vA@Manifest{}) +> (fA@Manifest{}) (gA@Manifest{}) +> (flagA@Manifest{}) del_t = +> let f'@Manifest{} = force $ traverse (Data.Array.Repa.zipWith (:*:) flagA +> (Data.Array.Repa.zipWith (:*:) uA +> (Data.Array.Repa.zipWith (:*:) vA fA))) id updatef +> flag = fsta +> u = fsta . snda +> v = fsta . snda . snda +> comp = snda . snda. snda +> updatef get c@(sh :. j :. i) = + +> if ((inBounds i j) && (i<=(imax-1))) then +> if (((flag $ get c) .&. _cf /= 0) +> && ((flag $ get (sh :. j :. (i+1))) .&. _cf /= 0)) then +> let + +> du2dx = ((((u $ get (sh :. j :. i)) +> +(u $ get (sh :. j :. (i+1))))**2) +> +gamma*(abs ((u $ get (sh :. j :. i)) +> +(u $ get (sh :. j :. (i+1))))) +> *((u $ get (sh :. j :. i)) +> -(u $ get (sh :. j :. (i+1)))) +> -(((u $ get (sh :. j :. (i-1))) +> +(u $ get (sh :. j :. i)))**2) +> -gamma*(abs ((u $ get (sh :. j :. (i-1))) +> +(u $ get (sh :. j :. i)))) +> *((u $ get (sh :. j :. (i-1))) +> -(u $ get (sh :. j :. i))))/(4.0*delx) +> duvdy = ((((v $ get (sh :. j :. i)) +> +(v $ get (sh :. j :. (i+1)))) +> *((u $ get (sh :. j :. i)) +> +(u $ get (sh :. (j+1) :. i)))) +> +(gamma*(abs ((v $ get (sh :. j :. i)) +> +(v $ get (sh :. j :. (i+1))))) +> *((u $ get (sh :. j :. i)) +> -(u $ get (sh :. (j+1) :. i)))) +> -(((v $ get (sh :. (j-1) :. i)) +> +(v $ get (sh :. (j-1) :. (i+1)))) +> *((u $ get (sh :. (j-1) :. i)) +> +(u $ get (sh :. j :. i)))) +> -(gamma*(abs ((v $ get (sh :. (j-1) :. i)) +> +(v $ get (sh :. (j-1) :. (i+1))))) +> *((u $ get (sh :. (j-1) :. i)) +> -(u $ get (sh :. j :. i)))))/(4.0*dely) +> laplu = ((u $ get (sh :. j :. (i+1))) +> -2.0*(u $ get (sh :. j :. i)) +> +(u $ get (sh :. j :. (i-1))))/delx/delx +> +((u $ get (sh :. (j+1) :. i)) +> -2.0*(u $ get (sh :. j :. i)) +> + (u $ get (sh :. (j-1) :. i)))/dely/dely +> in +> (u $ get (sh :. j :. i)) + del_t*(laplu/re-du2dx-duvdy) +> else +> u $ get (sh :. j :. i) +> else +> comp $ get (sh :. j :. i) +> g'@Manifest{} = force $ traverse (Data.Array.Repa.zipWith (:*:) flagA +> (Data.Array.Repa.zipWith (:*:) uA +> (Data.Array.Repa.zipWith (:*:) vA gA))) id updateg +> updateg get c@(sh :. j :. i) = +> if (inBounds i j) && (j<=(jmax-1)) then +> if (((flag $ get c) .&. _cf /= 0) +> && ((flag $ get (sh :. (j+1) :. i)) .&. _cf /= 0)) then +> let +> duvdx = ((((u $ get c)+(u $ get (sh :. (j+1) :. i))) +> *((v $ get c)+(v $ get (sh :. j :. (i+1))))) +> +(gamma*(abs ((u $ get c)+(u $ get (sh :. (j+1) :. i)))) +> *((v $ get c)-(v $ get (sh :. j :. (i+1))))) +> -(((u $ get (sh :. j :. (i-1))) +> +(u $ get (sh :. (j+1) :. (i-1)))) +> *((v $ get (sh :. j :. (i-1)))+(v $ get c))) +> -(gamma*(abs ((u $ get (sh :. j :. (i-1))) +> +(u $ get (sh :. (j+1) :. (i-1))))) +> *((v $ get (sh :. j :. (i-1)))-(v $ get c)))) +> /(4.0*delx) +> dv2dy = ((((v $ get c)+(v $ get (sh :. (j+1) :. i)))**2) +> +gamma*(abs ((v $ get c)+(v $ get (sh :. (j+1) :. i)))) +> *((v $ get c)-(v $ get (sh :. (j+1) :. i))) +> -(((v $ get (sh :. (j-1) :. i))+(v $ get c))**2) +> -gamma*(abs ((v $ get (sh :. (j-1) :. i))+(v $ get c))) +> *((v $ get (sh :. (j-1) :. i))-(v $ get c)))/(4.0*dely) +> laplv = ((v $ get (sh :. j :. (i+1))) +> -2*(v $ get c) +> +(v $ get (sh :. j :. (i-1))))/delx/delx + +> ((v $ get (sh :. (j+1) :. i)) +> -2*(v $ get c) +> +(v $ get (sh :. (j-1) :. i)))/dely/dely +> in +> (v $ get c) + del_t*(laplv/re-duvdx-dv2dy) +> else +> v $ get c +> else +> comp $ get c + +> f''@Manifest{} = force $ traverse (Data.Array.Repa.zipWith (:*:) uA f') id update +> where +> update get c@(sh :. j :. i) = +> if (j>=1 && j<=jmax) then +> if (i==0) then fsta $ get (sh :. j :. 0) +> else if (i==imax) then fsta $ get (sh :. j :. imax) +> else snda $ get (sh :. j :. i) +> else +> snda $ get (sh :. j :. i) +> g''@Manifest{} = force $ traverse (Data.Array.Repa.zipWith (:*:) vA g') id update +> where +> update get c@(sh :. j :. i) = +> if (i>=1 && i<=imax) then +> if (j==0) then fsta $ get (sh :. 0 :. i) +> else if (j==jmax) then fsta $ get (sh :. jmax :. i) +> else snda $ get (sh :. j :. i) +> else +> snda $ get (sh :. j :. i) +> in +> (f'', g'') +> + +> compute_rhs (f@Manifest{}) (g@Manifest{}) (rhs@Manifest{}) +> (flag@Manifest{}) del_t = +> force $ traverse (Data.Array.Repa.zipWith (:*:) flag +> (Data.Array.Repa.zipWith (:*:) rhs +> (Data.Array.Repa.zipWith (:*:) f g))) id update +> where +> update get c@(sh :. j :. i) = +> if (inBounds i j && (((flag $ get c) .&. _cf) /= 0)) then +> (((f $ get (sh :. j :. i))-(f $ get (sh :. j :. (i-1))))/delx + +> ((g $ get (sh :. j :. i))-(g $ get (sh :. (j-1) :. i)))/dely)/del_t +> else +> rhs $ get (sh :. j :. i) +> where +> flag = fsta +> rhs = fsta . snda +> f = fsta . snda . snda +> g = snda . snda . snda + + +> poisson :: Array DIM2 Double -> Array DIM2 Double -> Array DIM2 Int -> Int -> +> Double -> (Array DIM2 Double, Double, Int) +> poisson (p@Manifest{}) (rhs@Manifest{}) (flag@Manifest{}) ifluid res = +> let +> rdx2 = 1.0/(delx*delx) +> rdy2 = 1.0/(dely*dely) +> beta_2 = -omega/(2.0*(rdx2+rdy2)) + +> -- Calculate sum of squares +> sumSquares = snda $ toScalar (fold total (0 :*: 0.0) +> (fold square (0 :*: 0.0) (force $ (ignoreBoundary +> (Data.Array.Repa.zipWith (:*:) flag p))))) +> total (_ :*: y) (_ :*: x) = (0 :*: (y + x)) +> square (_ :*: y) x = if ((fsta x) .&. _cf /= 0) then +> 0 :*: (y+((snda x)**2)) +> else +> 0 :*: y + +> p0 = sqrt (sumSquares/((fromIntegral ifluid)::Double)) +> p0' = if (p0 < 0.0001) then 1.0 else p0 + +> -- Partial computation of residual +> computeRes (p@Manifest{}) p0 = +> let +> res'@Manifest{} = force $ traverse (Data.Array.Repa.zipWith (:*:) flag +> (Data.Array.Repa.zipWith (:*:) p rhs)) id update +> res = sumAll $ res' +> update get c@(sh :. j :. i) = +> let +> flag = fsta +> p = fsta . snda +> rhs = snda . snda +> in +> if ((inBounds i j) && ((flag $ get c) .&. _cf /= 0)) then +> ((((obstacle flag get (sh :. j :. (i+1)))* +> ((p $ get (sh :. j :. (i+1)))-(p $ get c))) +> -((obstacle flag get (sh :. j :. (i-1)))* +> ((p $ get c)-(p $ get (sh :. j :. (i-1))))))*rdx2 +> +(((obstacle flag get (sh :. (j+1) :. i))* +> ((p $ get (sh :. (j+1) :. i))-(p $ get c))) +> -((obstacle flag get (sh :. (j-1) :. i))* +> ((p $ get c) - (p $ get (sh :. (j-1) :. i)))))*rdy2 +> - (rhs $ get c))**2 +> else +> 0.0 +> in +> (sqrt (res/((fromIntegral ifluid)::Double)))/p0 + +> -- Red/Black SOR-iteration (compute new p) +> iterop rb (p@Manifest{}) (rhs@Manifest{}) (flag@Manifest{}) = +> force $ traverse (Data.Array.Repa.zipWith (:*:) flag +> (Data.Array.Repa.zipWith (:*:) p rhs)) +> id (update rb) +> update rb get c@(sh :. j :. i) = +> let +> flag = fsta +> p = fsta . snda +> rhs = snda . snda +> -- CSE the central element access +> !_c = get $ c +> in +> if ((inBounds i j) && (((i+j) `mod` 2) == rb)) then +> if ((flag _c) == (_cf .|. _bnsew)) then +> -- five point star for interior fluid cells +> (1.0-omega)*(p _c) - +> (beta_2 * +> (((p $ get (sh :. j :. (i+1))) +> +(p $ get (sh :. j :. (i-1))))*rdx2 +> + ((p $ get (sh :. (j+1) :. i)) +> +(p $ get (sh :. (j-1) :. i)))*rdy2 +> - (rhs _c))) +> else +> -- modified star near boundary +> if ((flag _c) .&. _cf /= 0) then +> let +> beta_mod = -omega/((obstacle flag get (sh :. j :. (i+1))+ +> (obstacle flag get (sh :. j :. (i-1))))*rdx2 +> +((obstacle flag get (sh :. (j+1) :. i))+ +> (obstacle flag get (sh :. (j-1) :. i)))*rdy2) +> in +> (1.0-omega)*(p _c) - (beta_mod*( +> (((obstacle flag get (sh :. j :. (i+1))*(p $ get (sh :. j :. (i+1))))) +> +((obstacle flag get (sh :. j :. (i-1))*(p $ get (sh :. j :. (i-1))))))*rdx2 +> +(((obstacle flag get (sh :. (j+1) :. i)*(p $ get (sh :. (j+1) :. i)))) +> +((obstacle flag get (sh :. (j-1) :. i)*(p $ get (sh :. (j-1) :. i)))))*rdy2 +> - (rhs _c))) +> else +> p _c +> else +> p _c + +> ----- CSE version (actually slower) + + let + _l = get (sh :. j :. (i-1)) + _r = get (sh :. j :. (i+1)) + _t = get (sh :. (j-1) :. i) + _b = get (sh :. (j+1) :. i) + beta_mod = -omega/(((obs $ flag _r)+ + (obs $ flag _l))*rdx2 + +((obs $ flag _b)+ + (obs $ flag _t))*rdy2) + in + (1.0-omega)*(p $ get c) - (beta_mod*( + ( + (((obs $ flag _r)*p _r) + +((obs $ flag _l)*p _l))*rdx2 + +(((obs $ flag _b)*p _b) + +((obs $ flag _t)*p _t))*rdy2 + -(rhs $ get c)))) + else + p $ get c + else + p $ get c + +> -- Iterations +> sorIterate iter (p@Manifest{}) res = +> let +> -- red +> (p'@Manifest{}) = iterop 0 p rhs flag +> -- black +> (p''@Manifest{}) = iterop 1 p' rhs flag +> -- res +> res' = computeRes p'' p0' +> in +> if (iter>=itermax) then +> (p, res, iter) +> else +> if (res' (p'', res', iter) +> else +> sorIterate (iter+1) p'' res' +> +> in +> sorIterate 0 p res + + +> update_velocity (u@Manifest{}) (v@Manifest{}) +> (f@Manifest{}) (g@Manifest{}) +> (p@Manifest{}) (flag@Manifest{}) del_t = +> let +> u'@Manifest{} = force $ traverse (Data.Array.Repa.zipWith (:*:) flag +> (Data.Array.Repa.zipWith (:*:) p +> (Data.Array.Repa.zipWith (:*:) f u))) id update +> where +> update get c@(sh :. j :. i) = +> let +> flag = fsta +> p = fsta . snda +> f = fsta . snda . snda +> u = snda . snda . snda +> in +> if (inBounds i j && (i<=(imax-1)) +> && ((flag $ get (sh :. j :. i)) .&. _cf /= 0) +> && ((flag $ get (sh :. j :. (i+1))) .&. _cf /= 0)) then +> (f $ get c)-((p $ get (sh :. j :. (i+1))) +> -(p $ get (sh :. j :. i)))*del_t/delx +> else +> u $ get c +> v'@Manifest{} = force $ traverse (Data.Array.Repa.zipWith (:*:) flag +> (Data.Array.Repa.zipWith (:*:) p +> (Data.Array.Repa.zipWith (:*:) g v))) id update +> where +> update get c@(sh :. j :. i) = +> let +> flag = fsta +> p = fsta . snda +> g = fsta . snda . snda +> v = snda . snda . snda +> in +> if (inBounds i j && (j<=(jmax-1)) +> && ((flag $ get (sh :. j :. i)) .&. _cf /= 0) +> && ((flag $ get (sh :. (j+1) :. i)) .&. _cf /= 0)) then +> (g $ get c)-((p $ get (sh :. (j+1) :. i)) +> -(p $ get (sh :. j :. i)))*del_t/dely +> else +> v $ get c +> in +> (u', v') \ No newline at end of file diff --git a/haskell/repa-2022/stack.yaml b/haskell/repa-2022/stack.yaml new file mode 100644 index 0000000..9ece693 --- /dev/null +++ b/haskell/repa-2022/stack.yaml @@ -0,0 +1,5 @@ +resolver: lts-18.18 +packages: [.] + +extra-deps: +- repa-3.4.1.5@sha256:4d94845e66d668345130f75f5de8f7c69139b91495c8c4b5a825014e8985e2fd,3562 diff --git a/haskell/repa-2022/stack.yaml.lock b/haskell/repa-2022/stack.yaml.lock new file mode 100644 index 0000000..64b22a5 --- /dev/null +++ b/haskell/repa-2022/stack.yaml.lock @@ -0,0 +1,19 @@ +# This file was autogenerated by Stack. +# You should not edit this file by hand. +# For more information, please see the documentation at: +# https://docs.haskellstack.org/en/stable/lock_files + +packages: +- completed: + hackage: repa-3.4.1.5@sha256:4d94845e66d668345130f75f5de8f7c69139b91495c8c4b5a825014e8985e2fd,3562 + pantry-tree: + size: 2965 + sha256: 1bb5ad1decd581b0d44b0e05cbdef5eb7967cdd80da481a06c3aab73e16087b9 + original: + hackage: repa-3.4.1.5@sha256:4d94845e66d668345130f75f5de8f7c69139b91495c8c4b5a825014e8985e2fd,3562 +snapshots: +- completed: + size: 586296 + url: https://raw.githubusercontent.com/commercialhaskell/stackage-snapshots/master/lts/18/18.yaml + sha256: 63539429076b7ebbab6daa7656cfb079393bf644971156dc349d7c0453694ac2 + original: lts-18.18 From 985c6882daee7e9e20d263b9a767430af3a96d7b Mon Sep 17 00:00:00 2001 From: Ben Orchard Date: Mon, 14 Feb 2022 11:17:52 +0000 Subject: [PATCH 2/3] update --- haskell/repa-2022/repa2022.cabal | 1 + haskell/repa-2022/src/Common.hs | 16 + haskell/repa-2022/src/Navier/Boundary.hs | 129 +++++++ haskell/repa-2022/src/Navier/Boundary.lhs | 122 ------- haskell/repa-2022/src/Navier/Datadefs.hs | 10 - haskell/repa-2022/src/Navier/Init.hs | 12 +- haskell/repa-2022/src/Navier/Output.hs | 1 + haskell/repa-2022/src/Navier/Simulation.hs | 357 ++++++++++++++++++++ haskell/repa-2022/src/Navier/Simulation.lhs | 351 ------------------- 9 files changed, 512 insertions(+), 487 deletions(-) create mode 100644 haskell/repa-2022/src/Common.hs create mode 100644 haskell/repa-2022/src/Navier/Boundary.hs delete mode 100644 haskell/repa-2022/src/Navier/Boundary.lhs create mode 100644 haskell/repa-2022/src/Navier/Simulation.hs delete mode 100644 haskell/repa-2022/src/Navier/Simulation.lhs diff --git a/haskell/repa-2022/repa2022.cabal b/haskell/repa-2022/repa2022.cabal index 82ddb99..6f0a7b1 100644 --- a/haskell/repa-2022/repa2022.cabal +++ b/haskell/repa-2022/repa2022.cabal @@ -10,6 +10,7 @@ build-type: Simple library exposed-modules: + Common Navier Navier.Boundary Navier.Datadefs diff --git a/haskell/repa-2022/src/Common.hs b/haskell/repa-2022/src/Common.hs new file mode 100644 index 0000000..b38f87c --- /dev/null +++ b/haskell/repa-2022/src/Common.hs @@ -0,0 +1,16 @@ +module Common where + +import Data.Array.Repa + +-- 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 arr = index arr Z diff --git a/haskell/repa-2022/src/Navier/Boundary.hs b/haskell/repa-2022/src/Navier/Boundary.hs new file mode 100644 index 0000000..dbe7359 --- /dev/null +++ b/haskell/repa-2022/src/Navier/Boundary.hs @@ -0,0 +1,129 @@ +{-# LANGUAGE NoMonomorphismRestriction #-} +{-# LANGUAGE QuasiQuotes #-} + +module Navier.Boundary where + +import Prelude hiding ( traverse ) + +import Data.Bits +import Data.Array.Repa +import Data.Array.Repa.Eval ( suspendedComputeP ) + +import Navier.Datadefs +import Navier.Params + +import Text.Printf +import Debug.Trace + +apply_boundary_conditions + :: Array U DIM2 Double -> Array U DIM2 Double -> Array U DIM2 Int + -> Int -> Int -> Double -> Double + -> (Array U DIM2 Double, Array U DIM2 Double) + +apply_boundary_conditions u v flag imax jmax ui vi = + + let ui' = traverse u id update + where update get c@(sh :. j :. i) + -- Fluid freely flows in from the west + | (i==0) = get (sh :. j :. 1) + -- Fluid freely flows out to the east + | (i==imax) = get (sh :. j :. (imax-1)) + | otherwise = get c + + u' = traverse ui' id update + where update get c@(sh :. j :. i) + -- fluid flows freely in the horizontal direction + | (j==(jmax+1)) = get (sh :. jmax :. i) + | (j==0) = get (sh :. 1 :. i) + | otherwise = get c + + vi' = traverse v id update + where update get c@(sh :. j :. i) + -- Fluid freely flows in from the west + | (i==0) = get (sh :. j :. 1) + -- Fluid freely flows out to the east + | (i==(imax+1)) = get (sh :. j :. imax) + | otherwise = get c + + v' = traverse vi' id update + where update get c@(sh :. j :. i) + -- vertical velocity approaches 0 at the north and south + | (j==jmax) = 0.0 + | (j==0) = 0.0 + | otherwise = get c + -- Apply no-slip boundary conditions to cells that are adjacent to + -- internal obstacle cells. This forces the u and v velocity to + -- tend towards zero in these cells. + u'' = traverse2 u' flag (\x y -> x) update + where + update getU getFlag c@(sh :. j :. i) = + if (inBounds i j && (getFlag c .&. _bnsew /= 0)) then + case (getFlag c) of + [p|_bn|] -> -getU (sh :. (j+1) :. i) + [p|_be|] -> 0.0 + [p|_bne|] -> 0.0 + [p|_bse|] -> 0.0 + [p|_bnw|] -> -getU (sh :. (j+1) :. i) + [p|_bs|] -> -getU (sh :. (j-1) :. i) + [p|_bsw|] -> -getU (sh :. (j-1) :. i) + _ -> getU c + else + getU c + u''' = traverse2 u'' flag (\x y -> x) update + where + update getU getFlag c@(sh :. j :. i) = + if (((inBounds i j || i==0) && (i<=(imax-1))) && + (getFlag (sh :. j :. (i+1)) .&. _bnsew /= 0)) then + case (getFlag (sh :. j :. (i+1))) of + [p|_bn|] -> -getU (sh :. (j+1) :. i) + [p|_bw|] -> 0.0 + [p|_bne|] -> -getU (sh :. (j+1) :. i) + [p|_bsw|] -> 0.0 + [p|_bnw|] -> 0.0 + [p|_bs|] -> -getU (sh :. (j-1) :. i) + [p|_bse|] -> -getU (sh :. (j-1) :. i) + _ -> getU c + else + getU c + v'' = traverse2 v' flag (\x y -> x) update + where + update getV getFlag c@(sh :. j :. i) = + if (inBounds i j && (getFlag c .&. _bnsew /= 0)) then + case (getFlag c) of + [p|_bn|] -> 0.0 + [p|_be|] -> -getV (sh :. j :. (i+1)) + [p|_bne|] -> 0.0 + [p|_bse|] -> -getV (sh :. j :. (i+1)) + [p|_bnw|] -> 0.0 + [p|_bw|] -> -getV (sh :. j :. (i-1)) + [p|_bsw|] -> -getV (sh :. j :. (i-1)) + _ -> getV c + else + getV c + v''' = traverse2 v'' flag (\x y -> x) update + where + update getV getFlag c@(sh :. j :. i) = + if (((inBounds i j || j==0) && j<=(jmax-1)) && + (getFlag (sh :. (j+1) :. i) .&. _bnsew /= 0)) then + case (getFlag (sh :. (j+1) :. i)) of + [p|_be|] -> -getV (sh :. j :. (i+1)) + [p|_bs|] -> 0.0 + [p|_bne|] -> -getV (sh :. j :. (i+1)) + [p|_bse|] -> 0.0 + [p|_bsw|] -> 0.0 + [p|_bw|] -> -getV (sh :. j :. (i-1)) + [p|_bnw|] -> -getV (sh :. j :. (i-1)) + _ -> getV c + else + getV c + -- Finally, fix the horizontal velocity at the western edge to have + -- a continual flow of fluid into the simulation. + u4 = suspendedComputeP $ traverse u''' id update + where update get c@(sh :. j :. i) + | (i==0 && j>=1 && j<=jmax) = ui + | otherwise = get c + v4 = suspendedComputeP $ traverse v''' id update + where update get c@(sh :. j :. i) + | (i==0 && j<=jmax) = 2*vi-(get (sh :. j :. 1)) + | otherwise = get c + in (u4, v4) diff --git a/haskell/repa-2022/src/Navier/Boundary.lhs b/haskell/repa-2022/src/Navier/Boundary.lhs deleted file mode 100644 index f838d85..0000000 --- a/haskell/repa-2022/src/Navier/Boundary.lhs +++ /dev/null @@ -1,122 +0,0 @@ -> {-# LANGUAGE NoMonomorphismRestriction #-} -> {-# LANGUAGE QuasiQuotes #-} - -> module Navier.Boundary where - -> import Data.Bits -> import Data.Array.Repa - -> import Navier.Datadefs -> import Navier.Params - -> import Text.Printf -> import Debug.Trace - -> apply_boundary_conditions :: Array DIM2 Double -> Array DIM2 Double -> -> Array DIM2 Int -> Int -> Int -> -> Double -> Double -> -> (Array DIM2 Double, Array DIM2 Double) -> apply_boundary_conditions (u@Manifest{}) (v@Manifest{}) -> (flag@Manifest{}) imax jmax ui vi = -> let ui'@Manifest{} = force $ traverse u id update -> where update get c@(sh :. j :. i) -> -- Fluid freely flows in from the west -> | (i==0) = get (sh :. j :. 1) -> -- Fluid freely flows out to the east -> | (i==imax) = get (sh :. j :. (imax-1)) -> | otherwise = get c -> u'@Manifest{} = force $ traverse ui' id update -> where update get c@(sh :. j :. i) -> -- fluid flows freely in the horizontal direction -> | (j==(jmax+1)) = get (sh :. jmax :. i) -> | (j==0) = get (sh :. 1 :. i) -> | otherwise = get c -> vi'@Manifest{} = force $ traverse v id update -> where update get c@(sh :. j :. i) -> -- Fluid freely flows in from the west -> | (i==0) = get (sh :. j :. 1) -> -- Fluid freely flows out to the east -> | (i==(imax+1)) = get (sh :. j :. imax) -> | otherwise = get c -> v'@Manifest{} = force $ traverse vi' id update -> where update get c@(sh :. j :. i) -> -- vertical velocity approaches 0 at the north and south -> | (j==jmax) = 0.0 -> | (j==0) = 0.0 -> | otherwise = get c -> -- Apply no-slip boundary conditions to cells that are adjacent to -> -- internal obstacle cells. This forces the u and v velocity to -> -- tend towards zero in these cells. -> u''@Manifest{} = force $ traverse2 u' flag (\x y -> x) update -> where -> update getU getFlag c@(sh :. j :. i) = -> if (inBounds i j && (getFlag c .&. _bnsew /= 0)) then -> case (getFlag c) of -> [p|_bn|] -> -getU (sh :. (j+1) :. i) -> [p|_be|] -> 0.0 -> [p|_bne|] -> 0.0 -> [p|_bse|] -> 0.0 -> [p|_bnw|] -> -getU (sh :. (j+1) :. i) -> [p|_bs|] -> -getU (sh :. (j-1) :. i) -> [p|_bsw|] -> -getU (sh :. (j-1) :. i) -> _ -> getU c -> else -> getU c -> u'''@Manifest{} = force $ traverse2 u'' flag (\x y -> x) update -> where -> update getU getFlag c@(sh :. j :. i) = -> if (((inBounds i j || i==0) && (i<=(imax-1))) && -> (getFlag (sh :. j :. (i+1)) .&. _bnsew /= 0)) then -> case (getFlag (sh :. j :. (i+1))) of -> [p|_bn|] -> -getU (sh :. (j+1) :. i) -> [p|_bw|] -> 0.0 -> [p|_bne|] -> -getU (sh :. (j+1) :. i) -> [p|_bsw|] -> 0.0 -> [p|_bnw|] -> 0.0 -> [p|_bs|] -> -getU (sh :. (j-1) :. i) -> [p|_bse|] -> -getU (sh :. (j-1) :. i) -> _ -> getU c -> else -> getU c -> v''@Manifest{} = force $ traverse2 v' flag (\x y -> x) update -> where -> update getV getFlag c@(sh :. j :. i) = -> if (inBounds i j && (getFlag c .&. _bnsew /= 0)) then -> case (getFlag c) of -> [p|_bn|] -> 0.0 -> [p|_be|] -> -getV (sh :. j :. (i+1)) -> [p|_bne|] -> 0.0 -> [p|_bse|] -> -getV (sh :. j :. (i+1)) -> [p|_bnw|] -> 0.0 -> [p|_bw|] -> -getV (sh :. j :. (i-1)) -> [p|_bsw|] -> -getV (sh :. j :. (i-1)) -> _ -> getV c -> else -> getV c -> v'''@Manifest{} = force $ traverse2 v'' flag (\x y -> x) update -> where -> update getV getFlag c@(sh :. j :. i) = -> if (((inBounds i j || j==0) && j<=(jmax-1)) && -> (getFlag (sh :. (j+1) :. i) .&. _bnsew /= 0)) then -> case (getFlag (sh :. (j+1) :. i)) of -> [p|_be|] -> -getV (sh :. j :. (i+1)) -> [p|_bs|] -> 0.0 -> [p|_bne|] -> -getV (sh :. j :. (i+1)) -> [p|_bse|] -> 0.0 -> [p|_bsw|] -> 0.0 -> [p|_bw|] -> -getV (sh :. j :. (i-1)) -> [p|_bnw|] -> -getV (sh :. j :. (i-1)) -> _ -> getV c -> else -> getV c -> -- Finally, fix the horizontal velocity at the western edge to have -> -- a continual flow of fluid into the simulation. -> u4@Manifest{} = force $ traverse u''' id update -> where update get c@(sh :. j :. i) -> | (i==0 && j>=1 && j<=jmax) = ui -> | otherwise = get c -> v4@Manifest{} = force $ traverse v''' id update -> where update get c@(sh :. j :. i) -> | (i==0 && j<=jmax) = 2*vi-(get (sh :. j :. 1)) -> | otherwise = get c -> in (u4, v4) \ No newline at end of file diff --git a/haskell/repa-2022/src/Navier/Datadefs.hs b/haskell/repa-2022/src/Navier/Datadefs.hs index fd25596..fd411a5 100644 --- a/haskell/repa-2022/src/Navier/Datadefs.hs +++ b/haskell/repa-2022/src/Navier/Datadefs.hs @@ -12,16 +12,6 @@ import Data.Array.Repa import Debug.Trace --- 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 - _cb = 0x0000 _bn = 0x0001 _bs = 0x0002 diff --git a/haskell/repa-2022/src/Navier/Init.hs b/haskell/repa-2022/src/Navier/Init.hs index d156b7c..1497cc0 100644 --- a/haskell/repa-2022/src/Navier/Init.hs +++ b/haskell/repa-2022/src/Navier/Init.hs @@ -6,6 +6,7 @@ import Prelude hiding ( traverse ) import Data.Bits import Data.Array.Repa +import Data.Array.Repa.Eval ( suspendedComputeP ) import Navier.Datadefs import Navier.Params @@ -14,12 +15,12 @@ init_flag :: Int -> Int -> Double -> Double -> (Array U DIM2 Int, Int) init_flag imax jmax delx dely = let dims = Z :. (jmax+2) :. (imax+2) - flag = suspendedComputeP $ fromList dims (Prelude.replicate ((imax+2)*(jmax+2)) (0::Int)) + flag = fromListUnboxed dims (Prelude.replicate ((imax+2)*(jmax+2)) (0::Int)) -- Make a circular obstacle mx = 20.0/41.0*(intCast jmax)*dely my = mx rad1 = 5.0/41.0*(intCast jmax)*dely - flag' = suspendedComputeP $ traverse flag id update + flag' = traverse flag id update where update get c@(sh :. j :. i) = if (inBounds i j) then @@ -30,13 +31,13 @@ init_flag imax jmax delx dely = else get c -- Make the north, sourth, east, and west boundary cells - flagi' = suspendedComputeP $ traverse flag' id update + flagi' = traverse flag' id update where update get c@(sh :. j :. i) = if (j==0) then _cb else if (j==(jmax+1)) then _cb else get c - flag'' = suspendedComputeP $ traverse flagi' id update + flag'' = traverse flagi' id update where update get c@(sh :. j :. i) = if (j>=1 && j<=jmax) then @@ -69,3 +70,6 @@ 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 diff --git a/haskell/repa-2022/src/Navier/Output.hs b/haskell/repa-2022/src/Navier/Output.hs index 69e466f..0124aa3 100644 --- a/haskell/repa-2022/src/Navier/Output.hs +++ b/haskell/repa-2022/src/Navier/Output.hs @@ -11,6 +11,7 @@ import System.IO import Text.Printf import Navier.Datadefs +import Common import Navier.Params type ArrayU a = Array U DIM2 a diff --git a/haskell/repa-2022/src/Navier/Simulation.hs b/haskell/repa-2022/src/Navier/Simulation.hs new file mode 100644 index 0000000..4e47b0a --- /dev/null +++ b/haskell/repa-2022/src/Navier/Simulation.hs @@ -0,0 +1,357 @@ +{-# LANGUAGE QuasiQuotes #-} +{-# LANGUAGE TypeOperators #-} + +module Navier.Simulation where + +import Prelude hiding ( traverse ) + +import Navier.Datadefs +import Common +import Navier.Params +import Data.Array.Repa +import Data.Array.Repa.Eval ( suspendedComputeP ) + +import Debug.Trace +import Data.Bits + +import Text.Printf + +set_timestamp_interval del_t u v = + -- del_t satisfying CFL conditions + -- else no time stepsize control + if (tau >= 1.0e-10) then + let + umax = foldAllS (\x y -> max (abs x) (abs y)) (1.0e-10) + (backpermute (Z :. (jmax-1) :. imax) + (\(sh :. j :. i) -> (sh :. (j+1) :. i)) u) + vmax = foldAllS (\x y -> max (abs x) (abs y)) (1.0e-10) + (backpermute (Z :. jmax :. (imax-1)) + (\(sh :. j :. i) -> (sh :. j :. (i+1))) v) + deltu = delx/umax + deltv = dely/vmax + deltre = 1.0/(1.0/(delx*delx)+1/(dely*dely))*re/2.0 + in + tau * (if (deltu=1 && j<=jmax) then + if (i==0) then fsta $ get (sh :. j :. 0) + else if (i==imax) then fsta $ get (sh :. j :. imax) + else snda $ get (sh :. j :. i) + else + snda $ get (sh :. j :. i) + g'' = suspendedComputeP $ traverse (Data.Array.Repa.zipWith (:*:) vA g') id update + where + update get c@(sh :. j :. i) = + if (i>=1 && i<=imax) then + if (j==0) then fsta $ get (sh :. 0 :. i) + else if (j==jmax) then fsta $ get (sh :. jmax :. i) + else snda $ get (sh :. j :. i) + else + snda $ get (sh :. j :. i) + in + (f'', g'') + +compute_rhs f g rhs flag del_t = + suspendedComputeP $ traverse (Data.Array.Repa.zipWith (:*:) flag + (Data.Array.Repa.zipWith (:*:) rhs + (Data.Array.Repa.zipWith (:*:) f g))) id update + where + update get c@(sh :. j :. i) = + if (inBounds i j && (((flag $ get c) .&. _cf) /= 0)) then + (((f $ get (sh :. j :. i))-(f $ get (sh :. j :. (i-1))))/delx + + ((g $ get (sh :. j :. i))-(g $ get (sh :. (j-1) :. i)))/dely)/del_t + else + rhs $ get (sh :. j :. i) + where + flag = fsta + rhs = fsta . snda + f = fsta . snda . snda + g = snda . snda . snda + + +poisson + :: Array U DIM2 Double -> Array U DIM2 Double -> Array U DIM2 Int + -> Int -> Double + -> (Array U DIM2 Double, Double, Int) +poisson p rhs flag ifluid res = + let + rdx2 = 1.0/(delx*delx) + rdy2 = 1.0/(dely*dely) + beta_2 = -omega/(2.0*(rdx2+rdy2)) + + -- Calculate sum of squares + sumSquares = + snda $ + toScalar + (foldS total (0 :*: 0.0) + (foldS square (0 :*: 0.0) + (ignoreBoundary (Data.Array.Repa.zipWith (:*:) flag p)))) + total (_ :*: y) (_ :*: x) = (0 :*: (y + x)) + square (_ :*: y) x = if ((fsta x) .&. _cf /= 0) then + 0 :*: (y+((snda x)**2)) + else + 0 :*: y + + p0 = sqrt (sumSquares/((fromIntegral ifluid)::Double)) + p0' = if (p0 < 0.0001) then 1.0 else p0 + + -- Partial computation of residual + computeRes p p0 = + let + res' = traverse (Data.Array.Repa.zipWith (:*:) flag + (Data.Array.Repa.zipWith (:*:) p rhs)) id update + res = sumAllS $ res' + update get c@(sh :. j :. i) = + let + flag = fsta + p = fsta . snda + rhs = snda . snda + in + if ((inBounds i j) && ((flag $ get c) .&. _cf /= 0)) then + ((((obstacle flag get (sh :. j :. (i+1)))* + ((p $ get (sh :. j :. (i+1)))-(p $ get c))) + -((obstacle flag get (sh :. j :. (i-1)))* + ((p $ get c)-(p $ get (sh :. j :. (i-1))))))*rdx2 + +(((obstacle flag get (sh :. (j+1) :. i))* + ((p $ get (sh :. (j+1) :. i))-(p $ get c))) + -((obstacle flag get (sh :. (j-1) :. i))* + ((p $ get c) - (p $ get (sh :. (j-1) :. i)))))*rdy2 + - (rhs $ get c))**2 + else + 0.0 + in + (sqrt (res/((fromIntegral ifluid)::Double)))/p0 + + -- Red/Black SOR-iteration (compute new p) + iterop rb p rhs flag = + traverse (Data.Array.Repa.zipWith (:*:) flag + (Data.Array.Repa.zipWith (:*:) p rhs)) + id (update rb) + update rb get c@(sh :. j :. i) = + let + flag = fsta + p = fsta . snda + rhs = snda . snda + -- CSE the central element access + !_c = get $ c + in + if ((inBounds i j) && (((i+j) `mod` 2) == rb)) then + if ((flag _c) == (_cf .|. _bnsew)) then + -- five point star for interior fluid cells + (1.0-omega)*(p _c) - + (beta_2 * + (((p $ get (sh :. j :. (i+1))) + +(p $ get (sh :. j :. (i-1))))*rdx2 + + ((p $ get (sh :. (j+1) :. i)) + +(p $ get (sh :. (j-1) :. i)))*rdy2 + - (rhs _c))) + else + -- modified star near boundary + if ((flag _c) .&. _cf /= 0) then + let + beta_mod = -omega/((obstacle flag get (sh :. j :. (i+1))+ + (obstacle flag get (sh :. j :. (i-1))))*rdx2 + +((obstacle flag get (sh :. (j+1) :. i))+ + (obstacle flag get (sh :. (j-1) :. i)))*rdy2) + in + (1.0-omega)*(p _c) - (beta_mod*( + (((obstacle flag get (sh :. j :. (i+1))*(p $ get (sh :. j :. (i+1))))) + +((obstacle flag get (sh :. j :. (i-1))*(p $ get (sh :. j :. (i-1))))))*rdx2 + +(((obstacle flag get (sh :. (j+1) :. i)*(p $ get (sh :. (j+1) :. i)))) + +((obstacle flag get (sh :. (j-1) :. i)*(p $ get (sh :. (j-1) :. i)))))*rdy2 + - (rhs _c))) + else + p _c + else + p _c + +{- + + ----- CSE version (actually slower) + + let + _l = get (sh :. j :. (i-1)) + _r = get (sh :. j :. (i+1)) + _t = get (sh :. (j-1) :. i) + _b = get (sh :. (j+1) :. i) + beta_mod = -omega/(((obs $ flag _r)+ + (obs $ flag _l))*rdx2 + +((obs $ flag _b)+ + (obs $ flag _t))*rdy2) + in + (1.0-omega)*(p $ get c) - (beta_mod*( + ( + (((obs $ flag _r)*p _r) + +((obs $ flag _l)*p _l))*rdx2 + +(((obs $ flag _b)*p _b) + +((obs $ flag _t)*p _t))*rdy2 + -(rhs $ get c)))) + else + p $ get c + else + p $ get c + +-} + + -- Iterations + sorIterate :: Int -> Array U DIM2 Double -> Double -> (Array U DIM2 Double, Double, Int) + sorIterate iter p res = + let + -- red + p' = iterop 0 p rhs flag + -- black + p'' = iterop 1 p' rhs flag + -- res + res' = computeRes p'' p0' + in + if (iter>=itermax) then + (p, res, iter) + else + if (res' module Navier.Simulation where - -> import Navier.Datadefs -> import Navier.Params -> import Data.Array.Repa - -> import Debug.Trace -> import Data.Bits - -> import Text.Printf - -> import Data.Array.Parallel.Base ((:*:)(..)) - -> set_timestamp_interval del_t (u@Manifest{}) (v@Manifest{}) = -> -- del_t satisfying CFL conditions -> -- else no time stepsize control -> if (tau >= 1.0e-10) then -> let -> umax = foldAll (\x y -> max (abs x) (abs y)) (1.0e-10) -> (backpermute (Z :. (jmax-1) :. imax) -> (\(sh :. j :. i) -> (sh :. (j+1) :. i)) u) -> vmax = foldAll (\x y -> max (abs x) (abs y)) (1.0e-10) -> (backpermute (Z :. jmax :. (imax-1)) -> (\(sh :. j :. i) -> (sh :. j :. (i+1))) v) -> deltu = delx/umax -> deltv = dely/vmax -> deltre = 1.0/(1.0/(delx*delx)+1/(dely*dely))*re/2.0 -> in -> tau * (if (deltu else min deltv deltre) -> else -> del_t - -> compute_tentative_velocity (uA@Manifest{}) (vA@Manifest{}) -> (fA@Manifest{}) (gA@Manifest{}) -> (flagA@Manifest{}) del_t = -> let f'@Manifest{} = force $ traverse (Data.Array.Repa.zipWith (:*:) flagA -> (Data.Array.Repa.zipWith (:*:) uA -> (Data.Array.Repa.zipWith (:*:) vA fA))) id updatef -> flag = fsta -> u = fsta . snda -> v = fsta . snda . snda -> comp = snda . snda. snda -> updatef get c@(sh :. j :. i) = - -> if ((inBounds i j) && (i<=(imax-1))) then -> if (((flag $ get c) .&. _cf /= 0) -> && ((flag $ get (sh :. j :. (i+1))) .&. _cf /= 0)) then -> let - -> du2dx = ((((u $ get (sh :. j :. i)) -> +(u $ get (sh :. j :. (i+1))))**2) -> +gamma*(abs ((u $ get (sh :. j :. i)) -> +(u $ get (sh :. j :. (i+1))))) -> *((u $ get (sh :. j :. i)) -> -(u $ get (sh :. j :. (i+1)))) -> -(((u $ get (sh :. j :. (i-1))) -> +(u $ get (sh :. j :. i)))**2) -> -gamma*(abs ((u $ get (sh :. j :. (i-1))) -> +(u $ get (sh :. j :. i)))) -> *((u $ get (sh :. j :. (i-1))) -> -(u $ get (sh :. j :. i))))/(4.0*delx) -> duvdy = ((((v $ get (sh :. j :. i)) -> +(v $ get (sh :. j :. (i+1)))) -> *((u $ get (sh :. j :. i)) -> +(u $ get (sh :. (j+1) :. i)))) -> +(gamma*(abs ((v $ get (sh :. j :. i)) -> +(v $ get (sh :. j :. (i+1))))) -> *((u $ get (sh :. j :. i)) -> -(u $ get (sh :. (j+1) :. i)))) -> -(((v $ get (sh :. (j-1) :. i)) -> +(v $ get (sh :. (j-1) :. (i+1)))) -> *((u $ get (sh :. (j-1) :. i)) -> +(u $ get (sh :. j :. i)))) -> -(gamma*(abs ((v $ get (sh :. (j-1) :. i)) -> +(v $ get (sh :. (j-1) :. (i+1))))) -> *((u $ get (sh :. (j-1) :. i)) -> -(u $ get (sh :. j :. i)))))/(4.0*dely) -> laplu = ((u $ get (sh :. j :. (i+1))) -> -2.0*(u $ get (sh :. j :. i)) -> +(u $ get (sh :. j :. (i-1))))/delx/delx -> +((u $ get (sh :. (j+1) :. i)) -> -2.0*(u $ get (sh :. j :. i)) -> + (u $ get (sh :. (j-1) :. i)))/dely/dely -> in -> (u $ get (sh :. j :. i)) + del_t*(laplu/re-du2dx-duvdy) -> else -> u $ get (sh :. j :. i) -> else -> comp $ get (sh :. j :. i) -> g'@Manifest{} = force $ traverse (Data.Array.Repa.zipWith (:*:) flagA -> (Data.Array.Repa.zipWith (:*:) uA -> (Data.Array.Repa.zipWith (:*:) vA gA))) id updateg -> updateg get c@(sh :. j :. i) = -> if (inBounds i j) && (j<=(jmax-1)) then -> if (((flag $ get c) .&. _cf /= 0) -> && ((flag $ get (sh :. (j+1) :. i)) .&. _cf /= 0)) then -> let -> duvdx = ((((u $ get c)+(u $ get (sh :. (j+1) :. i))) -> *((v $ get c)+(v $ get (sh :. j :. (i+1))))) -> +(gamma*(abs ((u $ get c)+(u $ get (sh :. (j+1) :. i)))) -> *((v $ get c)-(v $ get (sh :. j :. (i+1))))) -> -(((u $ get (sh :. j :. (i-1))) -> +(u $ get (sh :. (j+1) :. (i-1)))) -> *((v $ get (sh :. j :. (i-1)))+(v $ get c))) -> -(gamma*(abs ((u $ get (sh :. j :. (i-1))) -> +(u $ get (sh :. (j+1) :. (i-1))))) -> *((v $ get (sh :. j :. (i-1)))-(v $ get c)))) -> /(4.0*delx) -> dv2dy = ((((v $ get c)+(v $ get (sh :. (j+1) :. i)))**2) -> +gamma*(abs ((v $ get c)+(v $ get (sh :. (j+1) :. i)))) -> *((v $ get c)-(v $ get (sh :. (j+1) :. i))) -> -(((v $ get (sh :. (j-1) :. i))+(v $ get c))**2) -> -gamma*(abs ((v $ get (sh :. (j-1) :. i))+(v $ get c))) -> *((v $ get (sh :. (j-1) :. i))-(v $ get c)))/(4.0*dely) -> laplv = ((v $ get (sh :. j :. (i+1))) -> -2*(v $ get c) -> +(v $ get (sh :. j :. (i-1))))/delx/delx + -> ((v $ get (sh :. (j+1) :. i)) -> -2*(v $ get c) -> +(v $ get (sh :. (j-1) :. i)))/dely/dely -> in -> (v $ get c) + del_t*(laplv/re-duvdx-dv2dy) -> else -> v $ get c -> else -> comp $ get c - -> f''@Manifest{} = force $ traverse (Data.Array.Repa.zipWith (:*:) uA f') id update -> where -> update get c@(sh :. j :. i) = -> if (j>=1 && j<=jmax) then -> if (i==0) then fsta $ get (sh :. j :. 0) -> else if (i==imax) then fsta $ get (sh :. j :. imax) -> else snda $ get (sh :. j :. i) -> else -> snda $ get (sh :. j :. i) -> g''@Manifest{} = force $ traverse (Data.Array.Repa.zipWith (:*:) vA g') id update -> where -> update get c@(sh :. j :. i) = -> if (i>=1 && i<=imax) then -> if (j==0) then fsta $ get (sh :. 0 :. i) -> else if (j==jmax) then fsta $ get (sh :. jmax :. i) -> else snda $ get (sh :. j :. i) -> else -> snda $ get (sh :. j :. i) -> in -> (f'', g'') -> - -> compute_rhs (f@Manifest{}) (g@Manifest{}) (rhs@Manifest{}) -> (flag@Manifest{}) del_t = -> force $ traverse (Data.Array.Repa.zipWith (:*:) flag -> (Data.Array.Repa.zipWith (:*:) rhs -> (Data.Array.Repa.zipWith (:*:) f g))) id update -> where -> update get c@(sh :. j :. i) = -> if (inBounds i j && (((flag $ get c) .&. _cf) /= 0)) then -> (((f $ get (sh :. j :. i))-(f $ get (sh :. j :. (i-1))))/delx + -> ((g $ get (sh :. j :. i))-(g $ get (sh :. (j-1) :. i)))/dely)/del_t -> else -> rhs $ get (sh :. j :. i) -> where -> flag = fsta -> rhs = fsta . snda -> f = fsta . snda . snda -> g = snda . snda . snda - - -> poisson :: Array DIM2 Double -> Array DIM2 Double -> Array DIM2 Int -> Int -> -> Double -> (Array DIM2 Double, Double, Int) -> poisson (p@Manifest{}) (rhs@Manifest{}) (flag@Manifest{}) ifluid res = -> let -> rdx2 = 1.0/(delx*delx) -> rdy2 = 1.0/(dely*dely) -> beta_2 = -omega/(2.0*(rdx2+rdy2)) - -> -- Calculate sum of squares -> sumSquares = snda $ toScalar (fold total (0 :*: 0.0) -> (fold square (0 :*: 0.0) (force $ (ignoreBoundary -> (Data.Array.Repa.zipWith (:*:) flag p))))) -> total (_ :*: y) (_ :*: x) = (0 :*: (y + x)) -> square (_ :*: y) x = if ((fsta x) .&. _cf /= 0) then -> 0 :*: (y+((snda x)**2)) -> else -> 0 :*: y - -> p0 = sqrt (sumSquares/((fromIntegral ifluid)::Double)) -> p0' = if (p0 < 0.0001) then 1.0 else p0 - -> -- Partial computation of residual -> computeRes (p@Manifest{}) p0 = -> let -> res'@Manifest{} = force $ traverse (Data.Array.Repa.zipWith (:*:) flag -> (Data.Array.Repa.zipWith (:*:) p rhs)) id update -> res = sumAll $ res' -> update get c@(sh :. j :. i) = -> let -> flag = fsta -> p = fsta . snda -> rhs = snda . snda -> in -> if ((inBounds i j) && ((flag $ get c) .&. _cf /= 0)) then -> ((((obstacle flag get (sh :. j :. (i+1)))* -> ((p $ get (sh :. j :. (i+1)))-(p $ get c))) -> -((obstacle flag get (sh :. j :. (i-1)))* -> ((p $ get c)-(p $ get (sh :. j :. (i-1))))))*rdx2 -> +(((obstacle flag get (sh :. (j+1) :. i))* -> ((p $ get (sh :. (j+1) :. i))-(p $ get c))) -> -((obstacle flag get (sh :. (j-1) :. i))* -> ((p $ get c) - (p $ get (sh :. (j-1) :. i)))))*rdy2 -> - (rhs $ get c))**2 -> else -> 0.0 -> in -> (sqrt (res/((fromIntegral ifluid)::Double)))/p0 - -> -- Red/Black SOR-iteration (compute new p) -> iterop rb (p@Manifest{}) (rhs@Manifest{}) (flag@Manifest{}) = -> force $ traverse (Data.Array.Repa.zipWith (:*:) flag -> (Data.Array.Repa.zipWith (:*:) p rhs)) -> id (update rb) -> update rb get c@(sh :. j :. i) = -> let -> flag = fsta -> p = fsta . snda -> rhs = snda . snda -> -- CSE the central element access -> !_c = get $ c -> in -> if ((inBounds i j) && (((i+j) `mod` 2) == rb)) then -> if ((flag _c) == (_cf .|. _bnsew)) then -> -- five point star for interior fluid cells -> (1.0-omega)*(p _c) - -> (beta_2 * -> (((p $ get (sh :. j :. (i+1))) -> +(p $ get (sh :. j :. (i-1))))*rdx2 -> + ((p $ get (sh :. (j+1) :. i)) -> +(p $ get (sh :. (j-1) :. i)))*rdy2 -> - (rhs _c))) -> else -> -- modified star near boundary -> if ((flag _c) .&. _cf /= 0) then -> let -> beta_mod = -omega/((obstacle flag get (sh :. j :. (i+1))+ -> (obstacle flag get (sh :. j :. (i-1))))*rdx2 -> +((obstacle flag get (sh :. (j+1) :. i))+ -> (obstacle flag get (sh :. (j-1) :. i)))*rdy2) -> in -> (1.0-omega)*(p _c) - (beta_mod*( -> (((obstacle flag get (sh :. j :. (i+1))*(p $ get (sh :. j :. (i+1))))) -> +((obstacle flag get (sh :. j :. (i-1))*(p $ get (sh :. j :. (i-1))))))*rdx2 -> +(((obstacle flag get (sh :. (j+1) :. i)*(p $ get (sh :. (j+1) :. i)))) -> +((obstacle flag get (sh :. (j-1) :. i)*(p $ get (sh :. (j-1) :. i)))))*rdy2 -> - (rhs _c))) -> else -> p _c -> else -> p _c - -> ----- CSE version (actually slower) - - let - _l = get (sh :. j :. (i-1)) - _r = get (sh :. j :. (i+1)) - _t = get (sh :. (j-1) :. i) - _b = get (sh :. (j+1) :. i) - beta_mod = -omega/(((obs $ flag _r)+ - (obs $ flag _l))*rdx2 - +((obs $ flag _b)+ - (obs $ flag _t))*rdy2) - in - (1.0-omega)*(p $ get c) - (beta_mod*( - ( - (((obs $ flag _r)*p _r) - +((obs $ flag _l)*p _l))*rdx2 - +(((obs $ flag _b)*p _b) - +((obs $ flag _t)*p _t))*rdy2 - -(rhs $ get c)))) - else - p $ get c - else - p $ get c - -> -- Iterations -> sorIterate iter (p@Manifest{}) res = -> let -> -- red -> (p'@Manifest{}) = iterop 0 p rhs flag -> -- black -> (p''@Manifest{}) = iterop 1 p' rhs flag -> -- res -> res' = computeRes p'' p0' -> in -> if (iter>=itermax) then -> (p, res, iter) -> else -> if (res' (p'', res', iter) -> else -> sorIterate (iter+1) p'' res' -> -> in -> sorIterate 0 p res - - -> update_velocity (u@Manifest{}) (v@Manifest{}) -> (f@Manifest{}) (g@Manifest{}) -> (p@Manifest{}) (flag@Manifest{}) del_t = -> let -> u'@Manifest{} = force $ traverse (Data.Array.Repa.zipWith (:*:) flag -> (Data.Array.Repa.zipWith (:*:) p -> (Data.Array.Repa.zipWith (:*:) f u))) id update -> where -> update get c@(sh :. j :. i) = -> let -> flag = fsta -> p = fsta . snda -> f = fsta . snda . snda -> u = snda . snda . snda -> in -> if (inBounds i j && (i<=(imax-1)) -> && ((flag $ get (sh :. j :. i)) .&. _cf /= 0) -> && ((flag $ get (sh :. j :. (i+1))) .&. _cf /= 0)) then -> (f $ get c)-((p $ get (sh :. j :. (i+1))) -> -(p $ get (sh :. j :. i)))*del_t/delx -> else -> u $ get c -> v'@Manifest{} = force $ traverse (Data.Array.Repa.zipWith (:*:) flag -> (Data.Array.Repa.zipWith (:*:) p -> (Data.Array.Repa.zipWith (:*:) g v))) id update -> where -> update get c@(sh :. j :. i) = -> let -> flag = fsta -> p = fsta . snda -> g = fsta . snda . snda -> v = snda . snda . snda -> in -> if (inBounds i j && (j<=(jmax-1)) -> && ((flag $ get (sh :. j :. i)) .&. _cf /= 0) -> && ((flag $ get (sh :. (j+1) :. i)) .&. _cf /= 0)) then -> (g $ get c)-((p $ get (sh :. (j+1) :. i)) -> -(p $ get (sh :. j :. i)))*del_t/dely -> else -> v $ get c -> in -> (u', v') \ No newline at end of file From 83cbc1cd18292753f1b7f0f9cb6fbcdc994c1fe6 Mon Sep 17 00:00:00 2001 From: Ben Orchard Date: Mon, 14 Feb 2022 13:23:54 +0000 Subject: [PATCH 3/3] works but SLOW lol --- haskell/repa-2022/app/Main.hs | 76 ++++++++++++ haskell/repa-2022/package.yaml | 32 +++-- haskell/repa-2022/repa2022.cabal | 49 +++++++- haskell/repa-2022/src/Common.hs | 13 +-- haskell/repa-2022/src/Navier.lhs | 72 ------------ haskell/repa-2022/src/Navier/Datadefs.hs | 2 +- haskell/repa-2022/src/Navier/Init.hs | 4 +- haskell/repa-2022/src/Navier/Output.hs | 10 +- haskell/repa-2022/src/Navier/Params.hs | 4 +- haskell/repa-2022/src/Navier/Simulation.hs | 130 +++++++++++---------- 10 files changed, 224 insertions(+), 168 deletions(-) create mode 100644 haskell/repa-2022/app/Main.hs delete mode 100644 haskell/repa-2022/src/Navier.lhs diff --git a/haskell/repa-2022/app/Main.hs b/haskell/repa-2022/app/Main.hs new file mode 100644 index 0000000..bdcc778 --- /dev/null +++ b/haskell/repa-2022/app/Main.hs @@ -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)] diff --git a/haskell/repa-2022/package.yaml b/haskell/repa-2022/package.yaml index a9c8c94..04bc82e 100644 --- a/haskell/repa-2022/package.yaml +++ b/haskell/repa-2022/package.yaml @@ -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 diff --git a/haskell/repa-2022/repa2022.cabal b/haskell/repa-2022/repa2022.cabal index 6f0a7b1..a56e12e 100644 --- a/haskell/repa-2022/repa2022.cabal +++ b/haskell/repa-2022/repa2022.cabal @@ -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 diff --git a/haskell/repa-2022/src/Common.hs b/haskell/repa-2022/src/Common.hs index b38f87c..1359acd 100644 --- a/haskell/repa-2022/src/Common.hs +++ b/haskell/repa-2022/src/Common.hs @@ -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 diff --git a/haskell/repa-2022/src/Navier.lhs b/haskell/repa-2022/src/Navier.lhs deleted file mode 100644 index 03fcf85..0000000 --- a/haskell/repa-2022/src/Navier.lhs +++ /dev/null @@ -1,72 +0,0 @@ -> {-# LANGUAGE NoMonomorphismRestriction #-} - -> module Navier 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 = fromList dims (Prelude.replicate ((imax+2)*(jmax+2)) ui) -> v = fromList dims (Prelude.replicate ((imax+2)*(jmax+2)) vi) -> f = fromList dims (Prelude.replicate ((imax+2)*(jmax+2)) 0.0) -> g = fromList dims (Prelude.replicate ((imax+2)*(jmax+2)) 0.0) -> rhs = fromList dims (Prelude.replicate ((imax+2)*(jmax+2)) 0.0) -> p = fromList 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 DIM2 Double -> Array DIM2 Double -> -> Array DIM2 Double -> Array DIM2 Double -> Array DIM2 Double -> Array DIM2 Double -> -> Array DIM2 Int -> Int -> Int -> Int -> -> [(Array DIM2 Double, Array DIM2 Double, Array 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)] diff --git a/haskell/repa-2022/src/Navier/Datadefs.hs b/haskell/repa-2022/src/Navier/Datadefs.hs index fd411a5..a289877 100644 --- a/haskell/repa-2022/src/Navier/Datadefs.hs +++ b/haskell/repa-2022/src/Navier/Datadefs.hs @@ -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)) diff --git a/haskell/repa-2022/src/Navier/Init.hs b/haskell/repa-2022/src/Navier/Init.hs index 1497cc0..c2a6dfe 100644 --- a/haskell/repa-2022/src/Navier/Init.hs +++ b/haskell/repa-2022/src/Navier/Init.hs @@ -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 diff --git a/haskell/repa-2022/src/Navier/Output.hs b/haskell/repa-2022/src/Navier/Output.hs index 0124aa3..cc4343d 100644 --- a/haskell/repa-2022/src/Navier/Output.hs +++ b/haskell/repa-2022/src/Navier/Output.hs @@ -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) && diff --git a/haskell/repa-2022/src/Navier/Params.hs b/haskell/repa-2022/src/Navier/Params.hs index 2a3c6ee..77cb45a 100644 --- a/haskell/repa-2022/src/Navier/Params.hs +++ b/haskell/repa-2022/src/Navier/Params.hs @@ -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 #-} diff --git a/haskell/repa-2022/src/Navier/Simulation.hs b/haskell/repa-2022/src/Navier/Simulation.hs index 4e47b0a..1948447 100644 --- a/haskell/repa-2022/src/Navier/Simulation.hs +++ b/haskell/repa-2022/src/Navier/Simulation.hs @@ -8,6 +8,7 @@ import Prelude hiding ( traverse ) import Navier.Datadefs import Common import Navier.Params +import qualified Data.Array.Repa as Repa import Data.Array.Repa import Data.Array.Repa.Eval ( suspendedComputeP ) @@ -37,13 +38,13 @@ set_timestamp_interval del_t u v = del_t compute_tentative_velocity uA vA fA gA flagA del_t = - let f' = traverse (Data.Array.Repa.zipWith (:*:) flagA - (Data.Array.Repa.zipWith (:*:) uA - (Data.Array.Repa.zipWith (:*:) vA fA))) id updatef - flag = fsta - u = fsta . snda - v = fsta . snda . snda - comp = snda . snda. snda + let f' = traverse (Data.Array.Repa.zipWith (,) flagA + (Data.Array.Repa.zipWith (,) uA + (Data.Array.Repa.zipWith (,) vA fA))) id updatef + flag = fst + u = fst . snd + v = fst . snd . snd + comp = snd . snd. snd updatef get c@(sh :. j :. i) = if ((inBounds i j) && (i<=(imax-1))) then @@ -91,9 +92,9 @@ compute_tentative_velocity uA vA fA gA flagA del_t = u $ get (sh :. j :. i) else comp $ get (sh :. j :. i) - g' = traverse (Data.Array.Repa.zipWith (:*:) flagA - (Data.Array.Repa.zipWith (:*:) uA - (Data.Array.Repa.zipWith (:*:) vA gA))) id updateg + g' = traverse (Data.Array.Repa.zipWith (,) flagA + (Data.Array.Repa.zipWith (,) uA + (Data.Array.Repa.zipWith (,) vA gA))) id updateg updateg get c@(sh :. j :. i) = if (inBounds i j) && (j<=(jmax-1)) then if (((flag $ get c) .&. _cf /= 0) @@ -129,31 +130,31 @@ compute_tentative_velocity uA vA fA gA flagA del_t = else comp $ get c - f'' = suspendedComputeP $ traverse (Data.Array.Repa.zipWith (:*:) uA f') id update + f'' = suspendedComputeP $ traverse (Data.Array.Repa.zipWith (,) uA f') id update where update get c@(sh :. j :. i) = if (j>=1 && j<=jmax) then - if (i==0) then fsta $ get (sh :. j :. 0) - else if (i==imax) then fsta $ get (sh :. j :. imax) - else snda $ get (sh :. j :. i) + if (i==0) then fst $ get (sh :. j :. 0) + else if (i==imax) then fst $ get (sh :. j :. imax) + else snd $ get (sh :. j :. i) else - snda $ get (sh :. j :. i) - g'' = suspendedComputeP $ traverse (Data.Array.Repa.zipWith (:*:) vA g') id update + snd $ get (sh :. j :. i) + g'' = suspendedComputeP $ traverse (Data.Array.Repa.zipWith (,) vA g') id update where update get c@(sh :. j :. i) = if (i>=1 && i<=imax) then - if (j==0) then fsta $ get (sh :. 0 :. i) - else if (j==jmax) then fsta $ get (sh :. jmax :. i) - else snda $ get (sh :. j :. i) + if (j==0) then fst $ get (sh :. 0 :. i) + else if (j==jmax) then fst $ get (sh :. jmax :. i) + else snd $ get (sh :. j :. i) else - snda $ get (sh :. j :. i) + snd $ get (sh :. j :. i) in (f'', g'') compute_rhs f g rhs flag del_t = - suspendedComputeP $ traverse (Data.Array.Repa.zipWith (:*:) flag - (Data.Array.Repa.zipWith (:*:) rhs - (Data.Array.Repa.zipWith (:*:) f g))) id update + suspendedComputeP $ traverse (Data.Array.Repa.zipWith (,) flag + (Data.Array.Repa.zipWith (,) rhs + (Data.Array.Repa.zipWith (,) f g))) id update where update get c@(sh :. j :. i) = if (inBounds i j && (((flag $ get c) .&. _cf) /= 0)) then @@ -162,10 +163,10 @@ compute_rhs f g rhs flag del_t = else rhs $ get (sh :. j :. i) where - flag = fsta - rhs = fsta . snda - f = fsta . snda . snda - g = snda . snda . snda + flag = fst + rhs = fst . snd + f = fst . snd . snd + g = snd . snd . snd poisson @@ -179,17 +180,19 @@ poisson p rhs flag ifluid res = beta_2 = -omega/(2.0*(rdx2+rdy2)) -- Calculate sum of squares + sumSquares :: Double sumSquares = - snda $ - toScalar - (foldS total (0 :*: 0.0) - (foldS square (0 :*: 0.0) - (ignoreBoundary (Data.Array.Repa.zipWith (:*:) flag p)))) - total (_ :*: y) (_ :*: x) = (0 :*: (y + x)) - square (_ :*: y) x = if ((fsta x) .&. _cf /= 0) then - 0 :*: (y+((snda x)**2)) - else - 0 :*: y + snd $ toScalar $ sumSquaresCompute sumSquaresInner + sumSquaresInner :: Array D DIM2 (Int, Double) + sumSquaresInner = ignoreBoundary $ Data.Array.Repa.zipWith (,) flag p + sumSquaresCompute :: Array D DIM2 (Int, Double) -> Array U DIM0 (Int, Double) + sumSquaresCompute = foldS total (0, 0.0) . foldS square (0, 0.0) + + total (_, y) (_, x) = (0, (y + x)) + square (_, y) x = + if ((fst x) .&. _cf /= 0) + then (0, (y+((snd x)**2))) + else (0, y) p0 = sqrt (sumSquares/((fromIntegral ifluid)::Double)) p0' = if (p0 < 0.0001) then 1.0 else p0 @@ -197,14 +200,14 @@ poisson p rhs flag ifluid res = -- Partial computation of residual computeRes p p0 = let - res' = traverse (Data.Array.Repa.zipWith (:*:) flag - (Data.Array.Repa.zipWith (:*:) p rhs)) id update + res' = traverse (Data.Array.Repa.zipWith (,) flag + (Data.Array.Repa.zipWith (,) p rhs)) id update res = sumAllS $ res' update get c@(sh :. j :. i) = let - flag = fsta - p = fsta . snda - rhs = snda . snda + flag = fst + p = fst . snd + rhs = snd . snd in if ((inBounds i j) && ((flag $ get c) .&. _cf /= 0)) then ((((obstacle flag get (sh :. j :. (i+1)))* @@ -223,14 +226,14 @@ poisson p rhs flag ifluid res = -- Red/Black SOR-iteration (compute new p) iterop rb p rhs flag = - traverse (Data.Array.Repa.zipWith (:*:) flag - (Data.Array.Repa.zipWith (:*:) p rhs)) + traverse (Data.Array.Repa.zipWith (,) flag + (Data.Array.Repa.zipWith (,) p rhs)) id (update rb) update rb get c@(sh :. j :. i) = let - flag = fsta - p = fsta . snda - rhs = snda . snda + flag = fst + p = fst . snd + rhs = snd . snd -- CSE the central element access !_c = get $ c in @@ -302,14 +305,15 @@ poisson p rhs flag ifluid res = p'' = iterop 1 p' rhs flag -- res res' = computeRes p'' p0' + p''' = suspendedComputeP p'' in if (iter>=itermax) then (p, res, iter) else if (res'