Skip to content
Snippets Groups Projects
Commit 66132e59 authored by martijn's avatar martijn
Browse files

Made a separate file for the constants with a nice variable declaration style... 

Made a separate file for the constants with a nice variable declaration style common amongst physicists. Also factored out the derived simulationn parameters and renamed Main to SaltMarsh
parent 2d9456ae
Branches dev-martijn
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
salt-marsh/src-v2/Constants.hs 0 → 100644
+ 221
− 0
View file @ 66132e59
module SaltMarsh.Constants where
--------------------------------------------------------------------------------
-- Hydrodynamic parameters --
--------------------------------------------------------------------------------
h0 = 0.02 -- 0.02
g = 9.81 -- 9.81
nn = 0.05 -- 5e-2
hCrit = 0.001 -- 1e-3
slope = 0.00 -- 2e-3
difU = 5e-2 -- 1e-2
--hIn = 0.00003 -- 1e-4/5e-5
{- Symbol Dimension Description
- ----------------------------------------------------------------------------
- h0 m Imposed initial water layer thickness (equals
- the homogeneous equilibrium state)
-
- g m s^-2 Gravitational acceleration constant
- nn s m^(-1/3) Gauckler-Manning friction coefficient
- (Chow, 1959);
- (Mariotti and Fagherazzi, 2012 JGR) => n=0.016
-
- hCrit m Critical water layer thickness (is always
- retained) in wetting-drying algorithm
-
- slope m m^-1 Subsediment plane bed slope
- difU m^2 s^-1 Turbulent eddy viscosity
- hIn m Rain (G:later redefined?)
-}
--------------------------------------------------------------------------------
-- Vegetation roughness parameters --
--------------------------------------------------------------------------------
cb = 20.0
cd = 2.0
--hv = 0.5
kv = 0.41
-- difD = 1e-4
{- Symbol Dimension Description
- ----------------------------------------------------------------------------
- cb Chezy roughness of the bed
- cd Chezy friction coefficient for maximal
- vegetation density
-
- hv Shoot length
- kv Von Karman constant
- difD Lateral expansion of vegetation (G: dfined
- later?)
-}
--------------------------------------------------------------------------------
-- Sedimentation parameters --
--------------------------------------------------------------------------------
s_in = 0.0010 -- 1e-3
e0 = 10.0 -- 3e-2
pE = 0.9 -- 0.5 / 0.9
--d0 = 1e-4 -- 5e-2
--pD = 1.0 -- [0..1]
{- Symbol Dimension Description
- ----------------------------------------------------------------------------
- s_in m s^-1 Sediment input
- e0 s^-1 "Background" erosion rate.
- pE - Fraction by which sediment erosion is reduced
- if algae are at carrying capacity
-
- d0 m^2 s^-1 Sediment diffusivity in absence of algae
- (G:defined later?)
-
- pD - Fraction by which sediment diffusivity is
- reduced if algae are at carrying capacity
-}
--------------------------------------------------------------------------------
-- Diatom growth and erosion parameters --
--------------------------------------------------------------------------------
rr = 0.060 -- 0.30
kk = 1.0 -- 1.0
qq = 0.5 -- 0.031
qs = 0.006 -- 0.3
ec = 10.0 -- 20.0
dc = 1e-3
{- Symbol Dimension Description
- ----------------------------------------------------------------------------
- rr s^-1 Growth rate algae
- kk g m^-2 Carrying capacity algae
- qq m Value of water level where algae loss is
- (approximately) half max (higher Qq = lower
- overall loss rate)
-
- qs - Fraction of initial (homogeneous) water layer
- thickness at which sedimentation rate is
- (approximately) half maximal
-
- ec - Conversion factor from sediment erosion to
- algae loss
-
- dc - Conversion factor from sediment to biomass
- diffusion
-}
--------------------------------------------------------------------------------
-- Gradients --
--------------------------------------------------------------------------------
-- Gradient 1 (D0)
gradient_D0 = False
d0 = 1e-4 -- 1e-1
d0min = 1e-4 -- 1e-1
d0max = 1e-4 -- 1e-1
{- Symbol Dimension Description
- ----------------------------------------------------------------------------
- d0 m^2 s^-1 In absence of algae, background sediment
- d0min m^2 s^-1 diffusivity varies spatially between d0min and
- d0max m^2 s^-1 d0max
-}
-- Gradient 2 (Hin)
gradient_Hin = False
hIn = 0.00010 -- 3e-5 [m] Rain
hInMin = 0.00001 -- 1e-4 [m] Rain
hInMax = 0.00010 -- 1e-4 [m] Rain
{- Symbol Dimension Description
- ----------------------------------------------------------------------------
- d0 m^2 s^-1 In absence of algae, background sediment
- d0min m^2 s^-1 diffusivity varies spatially between d0min and
- d0max m^2 s^-1 d0max
-}
-- Gradient 3 (DifD)
gradient_DifD = False
difD = 1e-4 -- 1e-4 # Lateral expansion of vegetation
difDmin = 1e-6 -- 1e-5
difDmax = 1e-1 -- 1e-2
{- Symbol Dimension Description
- ----------------------------------------------------------------------------
- difD m^2 s^-1 Lateral expansion of vegetation
- difDmin m^2 s^-1
- difDmax m^2 s^-1
-}
-- Gradient 4 (pD)
gradient_pD = False
pD = 1e0 -- 0.0 [-] Fraction by which sediment diffusivity is reduced when at carrying capacity
pDmin = 1e-1 :: Exp Float -- 1e-5
pDmax = 1e1 :: Exp Float -- 1e-2
pDcon = 2.0 -- Constant determines curvature of logarithmic pD-gradient
{- Symbol Dimension Description
- ----------------------------------------------------------------------------
- pD - Fraction by which sediment diffusivity is
- reduced when at carrying capacity
-
- pDmin -
- pDmax -
- pDcon - Constant determines curvature of logarithmic
- pD-gradient
-}
-- Gradient 5 (Hv)
gradient_Hv :: Bool
gradient_Hv = False
hv = 5e-1 -- Shoot length
hvmin = 1e-3 -- 1e-5
hvmax = 1e0 -- 1e-2
{- Symbol Dimension Description
- ----------------------------------------------------------------------------
- hv ? Shoot length
- hvmin ?
- hvmax ?
-}
--------------------------------------------------------------------------------
-- Simulation constants --
--------------------------------------------------------------------------------
endTime = 2000
numFrames = 100
dT = 0.01
lengthX, lengthY :: Exp Float
lengthX = 100.0
lengthY = 100.0
block_Size_X = 16
block_Size_Y = 16
block_Number_X = 32
block_Number_Y = 32
{- Symbol Default Value Description
- ----------------------------------------------------------------------------
- endTime 20 Total time
- numFrames 50 Number of times the figure is updated
- dT 0.0005 Time step size (frame interval time)
- lengthX 800 Size of the domain in physical dimensions
- lengthY 100 Size of the domain in physical dimensions
- block_Size_X 16 Thread block size(*)
- block_Size_Y 16 Thread block size
- block_Number_X 32 ?
- block_Number_Y 32 ?
-
- * We do not need it here, just to calculate matrix sizes
-}
salt-marsh/src-v2/Derived.hs 0 → 100644
+ 27
− 0
View file @ 66132e59
module SaltMarsh.Derived where
-- Module
import SaltMarsh.Constants
{- Derived simulation parameters
-
- XXX: this is not a good way to do it. The size of the grid should be
- taken directly from the 'shape' of the input array.
- * TLM 2019-04-02
-}
grid_Width = block_Size_X * block_Number_X -- Matrix A width
grid_Height = block_Size_Y * block_Number_Y -- Matrix A height
dX, dY :: Exp Float
dX = lengthX / grid_Width
dY = lengthY / grid_Height
{- Symbol Description
- ----------------------------------------------------------------------------
- grid_Width Matrix A width
- grid_Height Matrix A height
- dX Space step
- dY Space step
-}
salt-marsh/src-v2/SaltMarsh.hs 0 → 100644
+ 333
− 0
View file @ 66132e59
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE StandaloneDeriving #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE ViewPatterns #-}
{-# LANGUAGE NoMonomorphismRestriction #-}
{-
- Salt marsh simulation, code ported to Haskell/Accelerate from Johan van der
- Koppel's Python & OpenCL implementation
-}
module SaltMarsh where
-- Standard library
import Codec.Picture
import Prelude as P
import System.Console.AsciiProgress
import System.Environment
import Text.Printf
import Text.Read ( readMaybe )
-- Accelerate
import Data.Array.Accelerate
import Data.Array.Accelerate.LLVM.Native
import Data.Array.Accelerate.System.Random.MWC
import qualified Data.Array.Accelerate as A
-- Module
import SaltMarsh.Constants
import SaltMarsh.Derived
type DataPoint =
( Float -- u
, Float -- v
, Float -- h
, Float -- s
, Float -- d
)
--Precomputation
h_homo = h0
u_homo = sqrt slope * h_homo**(2.0 / 3.0) / nn -- Balance between downslope acceleration and friction
s_homo = s_in * u_homo * u_homo * (h0 - hCrit)/(qs + h0 - hCrit) / e0
difco :: Exp Float -> Exp Float -> Exp Float -> Exp Float
difco d1 d2 column = d0x * ds
where
dm = (d1 + d2) / 2
pDx :: Exp Float
pDx | gradient_pD = exp(( (log pDmax - log pDmin) * column / grid_Width + log pDmin) /
(grid_Width + log pDmin))
| otherwise = pD
ds :: Exp Float
ds = exp (-dm * pDx / kk)
d0x | gradient_D0 = exp( ((log d0max)-(log d0min))* column/grid_Width + log d0min)
| otherwise = d0
-- projections named so they correspond to the array names in the OpenCL code
u, v, s ,d,h :: (Exp DataPoint) -> Exp Float
u tupleExp = let (x,_,_,_,_) = (unlift tupleExp) :: (Exp Float,Exp Float,Exp Float,Exp Float,Exp Float) in x
v tupleExp = let (_,x,_,_,_) = (unlift tupleExp) :: (Exp Float,Exp Float,Exp Float,Exp Float,Exp Float) in x
h tupleExp = let (_,_,x,_,_) = (unlift tupleExp) :: (Exp Float,Exp Float,Exp Float,Exp Float,Exp Float) in x
s tupleExp = let (_,_,_,x,_) = (unlift tupleExp) :: (Exp Float,Exp Float,Exp Float,Exp Float,Exp Float) in x
d tupleExp = let (_,_,_,_,x) = (unlift tupleExp) :: (Exp Float,Exp Float,Exp Float,Exp Float,Exp Float) in x
-- 'b' never seems to change, so in contrast to the OpenCL code, it's not returned as result
-- column info bundled with b, as we only have stencil2 predefined, not stencil3
-- is it really necessary for the simulation to know the column number?
simulateUV :: Stencil3x3 (Float, Float) -> Stencil3x3 DataPoint -> Exp DataPoint
simulateUV ((_ , cbtop, _ ),
(cbleft, cbcurr, cbright),
(_, cbbot, _ ))
((_, top, _ ),
(left, curr, right ),
(_, bot, _ ))
= lift (u_new, v_new, h_new, s_new, d_new)
where
-- bundling the columns and b together is ugly - do we need to cols?
column = A.fst cbcurr
bleft = A.snd cbleft
bright = A.snd cbright
bbot = A.snd cbbot
btop = A.snd cbtop
hvx :: Exp Float
hvx | gradient_Hv = exp ((log hvmax - log hvmin) * A.fst cbcurr / grid_Width + log hvmin )
| otherwise = hv
uabs = sqrt(u curr * u curr + v curr * v curr)
h_curr = A.max (h curr) hCrit
ct :: Exp Float
ct = sqrt (1.0/(1.0/(cb * cb) + 1.0 / (2.0 * g) * cd * d curr * hvx)) +
sqrt g / kv * log (A.max (h_curr / hvx) 1.0)
du = -g * dO_dx
- u curr * d_dx u
- v curr * d_dy u
- g / (ct * ct) * uabs * u curr / h_curr
+ difU * d2_dxy2 u
dv = -g * dO_dy
- u curr * d_dx v
- v curr * d_dy v
- g / (ct * ct) * uabs * v curr / h_curr
+ difU * (d2_dxy2 v)
u_new = u curr + du * dT
v_new = v curr + dv * dT
d_dy pop = (pop bot - pop top)/ 2.0 / dY
d_dx z = (z right - z left)/ 2.0 / dX
dO_dx = ((h right + s right + bright) -
(h left + s left + bleft )) / 2.0 / dX
dO_dy = ((h bot + s bot + bbot) -
(h top + s top + btop) ) / 2.0 / dY
d2_dxy2 z = (z left + z right - 2.0 * z curr) / dX / dX +
(z top + z bot - 2.0 * z curr) / dY / dY
phi = 1.0
-- other state variables
dh = - d_uh_dx - d_vh_dy
drh | gradient_Hin = exp((log hInMax -log hInMin) *
(A.fst cbcurr) /grid_Width +
log hInMin)
| otherwise = hIn
h_eff = h_curr - hCrit
ds = s_in * h_eff / (qs + h_eff)
- e0 * (1 - pE * d curr/kk) * s curr * uabs * uabs * (g / (ct*ct))
+ (d2_dxy2_S s d)
difDx
| gradient_DifD = exp((log difDmax - log difDmin) *
A.fst cbcurr / grid_Width +
log difDmin)
| otherwise = difD
dD = rr * d curr * (1 - d curr / kk) * qq / (qq + h_eff) -
ec * d curr * uabs * uabs * (g/(ct * ct)) + difDx * d2_dxy2 d
h_new = h_curr + dT * (dh + drh)
s_new = s curr + dT * ds * phi
d_new = d curr + dT * dD * phi
d_uh_dx = (u right * h right - u left * h left) / 2.0 / dX
d_vh_dy = (v bot * h bot - v top * h top) / 2.0 / dY
d2_dxy2_S s d = difco (d right) (d curr) column / dX / dX * (s right - s curr) -
difco (d curr) (d left) column / dX / dX * (s curr - s left) +
difco (d bot) (d curr) column / dY / dY * (s bot - s curr) -
difco (d curr) (d top) column / dY / dY * (s curr - s top )
-- after each step, we adjust the boundaries (still not exactly the same result as
-- doing the boundaries properly straight away
simulationStep :: Acc (Array DIM2 DataPoint) -> Acc (Array DIM2 DataPoint)
simulationStep arr = fixBoundaries $ stencil2 simulateUV clamp b clamp arr
where
fixBoundaries :: Acc (Array DIM2 DataPoint) -> Acc (Array DIM2 DataPoint)
fixBoundaries arr = generate (shape arr) boundaryFn
where
(uarr, varr, harr, sarr, darr) = A.unzip5 arr
boundaryFn :: (Exp DIM2 -> Exp DataPoint)
boundaryFn i =
cond (r A.== 0)
(lift $ ( 2 * uarr A.! row1 - uarr A.! row2 :: Exp Float
, 2 * varr A.! row1 - varr A.! row2 :: Exp Float
, harr A.! neumannTop :: Exp Float
, 0 :: Exp Float
, darr A.! neumannTop :: Exp Float
)
)(
cond (r A.== grid_Height - 1)
(lift $ ( uarr A.! neumannBot
, - varr A.! neumannBot
, harr A.! neumannBot
, sarr A.! neumannBot
, darr A.! neumannBot
)
)(
cond (c A.== 0)
(lift $ (-uarr A.! neumannLeft
, varr A.! neumannLeft
, harr A.! neumannLeft
, sarr A.! neumannLeft
, darr A.! neumannLeft
)
)(
cond (c A.== grid_Width - 1)
(lift $ (-uarr A.! neumannRight
, varr A.! neumannRight
, harr A.! neumannRight
, sarr A.! neumannRight
, darr A.! neumannRight
)
)(
lift $ ( uarr A.! i
, varr A.! i
, harr A.! i
, sarr A.! i
, darr A.! i
)
))))
where
(_ :. r :. c) = unlift i :: (Exp Z :. Exp Int :. Exp Int)
row1 :: Exp DIM2
row1 = lift (Z :. (r+1) :. (0 :: Exp Int))
row2 :: Exp DIM2
row2 = lift (Z :. (r+2) :. (0 :: Exp Int))
neumannTop = lift (Z :. (1 :: Exp Int) :. c)
neumannBot = lift (Z :. (grid_Height - 2 :: Exp Int) :. c)
neumannLeft = lift (Z :. r :. (1 :: Exp Int))
neumannRight = lift (Z :. r :. (grid_Width - 2 :: Exp Int))
b :: Acc (Array DIM2 (Float, Float))
b = generate (lift $ Z :. (grid_Height :: Int) :.
(grid_Width :: Int)) bfn
where
bfn :: Exp DIM2 -> Exp (Float, Float)
bfn i = lift (A.fromIntegral c :: Exp Float, (1.0 - (A.fromIntegral c)/grid_Width) * slope *lengthX :: Exp Float)
where
(_ :. r :. c) = unlift i :: (Exp Z :. Exp Int :. Exp Int)
whileCalcFrame :: Acc (Array DIM0 Int, Array DIM2 DataPoint) -> Acc (Array DIM0 Int, Array DIM2 DataPoint)
whileCalcFrame initState = awhile cont next initState
where
cont :: Acc (Array DIM0 Int, Array DIM2 DataPoint) -> Acc (Scalar Bool)
cont st = unit ((afst st) ! (lift Z) A.> 0)
next :: Acc (Array DIM0 Int, Array DIM2 DataPoint) -> Acc (Array DIM0 Int, Array DIM2 DataPoint)
next state = lift (cnt', simulationStep arr)
where arr = asnd state
cnt' = A.map (+(-1)) (afst state)
main :: IO ()
main = do
-- args <- getArgs
-- let steps = if P.null args
-- then 10
-- else case (readMaybe $ head args) of
-- Just s -> s
-- _ -> 10
-- set up initial conditions and run loop
ds <- randomArray rand (Z :. grid_Height :. grid_Width)
let dp0 = runN dataPoints ds
go = runN whileCalcFrame
loop :: ProgressBar -> Int -> Matrix DataPoint -> IO ()
loop !pg !i !dp
| i P.>= numFrames = complete pg -- strictly speaking should not be necessary
| otherwise = do
let -- compute next frame
(_, !dp') = go (fromList Z [stepsPerFrame], dp)
-- let fl = P.map (\(_,_,_,s,_) -> s) $ A.toList dp'
-- writeFile "saltMarsh.txt" $ show fl
-- or, to write png:
writePng (printf "output/sediment-%04d.png" i) $ heatMap grid_Width grid_Height sedimentSelector dp'
writePng (printf "output/waterflow-%04d.png" i) $ heatMap grid_Width grid_Height waterflowSelector dp'
tick pg
loop pg (i+1) dp'
-- run simulation
displayConsoleRegions $ do
printf "Running saltMarsh simulation with %d frames of %d steps each\n" (numFrames::Int) stepsPerFrame
pg <- newProgressBar def { pgTotal = numFrames
, pgOnCompletion = Just "Complete after :elapsed seconds"
}
loop pg 0 dp0
where
stepsPerFrame :: Int
stepsPerFrame = P.floor $ endTime / (dT * numFrames)
dataPoints :: Acc (Array DIM2 Float) -> Acc (Array DIM2 DataPoint)
dataPoints = A.map f
where
f :: Exp Float -> Exp DataPoint
f x = lift ( u_homo :: Exp Float
, 0.0 :: Exp Float
, h_homo :: Exp Float
, s_homo :: Exp Float
, x
)
sedimentSelector = \(_,_,_,s,_) -> s
waterflowSelector = \(u,v,_,_,_) -> sqrt (u*u + v*v)
-- simple heat map assuming values between 0 and 0.45
heatMap :: Int -> Int -> (DataPoint -> Float) -> Array DIM2 DataPoint -> Image PixelRGB8
heatMap width height sel values = generateImage heatPixel 512 512
where
heatPixel x y = col $ val y x
val x y = P.min (P.floor $ (sel $ indexArray values (Z :. x :. y)) * 1133) 510
col v = PixelRGB8
(P.fromIntegral $ P.max 0 (v - 255))
(P.fromIntegral $ P.min v 255)
(P.fromIntegral $ 255 - (P.min 255 v))
rand :: DIM2 :~> Float
rand ix gen = do
v <- uniformR (0,1 :: Float) ix gen
return $ if v P.< 0.02 then 1 else 0
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment