diff --git a/gloss-examples.cabal b/gloss-examples.cabal --- a/gloss-examples.cabal +++ b/gloss-examples.cabal @@ -1,5 +1,5 @@ Name: gloss-examples -Version: 1.7.4.5 +Version: 1.7.5.1 License: MIT License-file: LICENSE Author: Ben Lippmeier @@ -21,7 +21,7 @@ Build-depends: base == 4.*, gloss == 1.7.*, - bytestring == 0.9.*, + bytestring >= 0.9 && < 0.11, bmp == 1.2.* Main-is: Main.hs hs-source-dirs: picture/Bitmap @@ -148,7 +148,7 @@ Build-depends: base == 4.*, gloss == 1.7.*, - containers >= 0.3 && <= 0.5, + containers >= 0.3 && <= 0.6, ghc-prim == 0.2.* Main-is: Main.hs other-modules: Actor Advance Collide Config Contact QuadTree World @@ -189,7 +189,7 @@ Build-depends: base == 4.*, gloss == 1.7.*, - gloss-raster == 1.7.4.* + gloss-raster == 1.7.5.* Main-is: Main.hs hs-source-dirs: raster/Crystal ghc-options: @@ -198,13 +198,16 @@ -funfolding-use-threshold1000 -funfolding-keeness-factor1000 -fllvm -optlo-O3 - + extensions: + PatternGuards + BangPatterns + Executable gloss-ray Build-depends: base == 4.*, gloss == 1.7.*, - gloss-raster == 1.7.4.* + gloss-raster == 1.7.5.* Main-is: Main.hs other-modules: Light Object Trace Vec3 World hs-source-dirs: raster/Ray @@ -220,7 +223,7 @@ Build-depends: base == 4.*, gloss == 1.7.*, - gloss-raster == 1.7.4.* + gloss-raster == 1.7.5.* Main-is: Main.hs hs-source-dirs: raster/Pulse ghc-options: @@ -235,7 +238,7 @@ Build-depends: base == 4.*, gloss == 1.7.*, - gloss-raster == 1.7.4.*, + gloss-raster == 1.7.5.*, vector == 0.9.*, ghc-prim Main-is: Main.hs @@ -245,7 +248,6 @@ -Odph -fno-liberate-case -funfolding-use-threshold1000 -funfolding-keeness-factor1000 - -fllvm -optlo-O3 Executable gloss-fluid @@ -293,7 +295,7 @@ Build-depends base == 4.*, gloss == 1.7.*, - repa >= 3.2.1 && < 3.3 + repa == 3.2.* Main-is: Main.hs other-modules: Solver hs-source-dirs: raster/Mandel diff --git a/raster/Crystal/Main.hs b/raster/Crystal/Main.hs --- a/raster/Crystal/Main.hs +++ b/raster/Crystal/Main.hs @@ -4,10 +4,106 @@ -- Based on code from: -- http://mainisusuallyafunction.blogspot.com/2011/10/quasicrystals-as-sums-of-waves-in-plane.html -- -{-# LANGUAGE BangPatterns #-} import Graphics.Gloss.Raster.Field import System.Environment +import System.Exit +import Data.Char +-- Main ----------------------------------------------------------------------- +main :: IO () +main + = do args <- getArgs + config <- parseArgs args defaultConfig + + let display + = case configFullScreen config of + True -> FullScreen (configSizeX config, configSizeY config) + False -> InWindow "Crystal" + (configSizeX config, configSizeY config) + (10, 10) + + let scale = fromIntegral $ configScale config + animateField display + (configZoom config, configZoom config) + (quasicrystal scale (configDegree config)) + + +-- Config --------------------------------------------------------------------- +data Config + = Config + { configSizeX :: Int + , configSizeY :: Int + , configFullScreen :: Bool + , configZoom :: Int + , configScale :: Int + , configDegree :: Int } + deriving Show + +defaultConfig :: Config +defaultConfig + = Config + { configSizeX = 800 + , configSizeY = 600 + , configFullScreen = False + , configZoom = 2 + , configScale = 30 + , configDegree = 5 } + + +parseArgs :: [String] -> Config -> IO Config +parseArgs args config + | [] <- args + = return config + + | "-fullscreen" : sizeX : sizeY : rest <- args + , all isDigit sizeX + , all isDigit sizeY + = parseArgs rest + $ config { configSizeX = read sizeX + , configSizeY = read sizeY + , configFullScreen = True } + + | "-window" : sizeX : sizeY : rest <- args + , all isDigit sizeX + , all isDigit sizeY + = parseArgs rest + $ config { configSizeX = read sizeX + , configSizeY = read sizeY + , configFullScreen = False } + + | "-zoom" : zoom : rest <- args + , all isDigit zoom + = parseArgs rest + $ config { configZoom = read zoom } + + | "-scale" : scale : rest <- args + , all isDigit scale + = parseArgs rest + $ config { configScale = read scale } + + | "-degree" : degree : rest <- args + , all isDigit degree + = parseArgs rest + $ config { configDegree = read degree } + + | otherwise + = do printUsage + exitWith $ ExitFailure 1 + + +printUsage :: IO () +printUsage + = putStr $ unlines + [ "quazicrystal [flags]" + , " -fullscreen sizeX sizeY Run full screen" + , " -window sizeX sizeY Run in a window (default 800, 600)" + , " -zoom Pixel replication factor (default 5)" + , " -scale Feature size of visualisation (default 30)" + , " -degree Number waves to sum for each point (default 5)" + , "" + , " You'll want to run this with +RTS -N to enable threads" ] + + -- Types ---------------------------------------------------------------------- -- | Angle in radians. type Angle = Float @@ -78,30 +174,4 @@ rampColor v = rawColor v (0.4 + (v * 0.6)) 1 1 - --- Main ----------------------------------------------------------------------- -main :: IO () -main - = do args <- getArgs - case args of - [] -> run 800 600 2 30 5 - - [sizeX, sizeY, zoom, scale, degree] - -> run (read sizeX) (read sizeY) (read zoom) (read scale) (read degree) - - _ -> putStr $ unlines - [ "quazicrystal " - , " sizeX, sizeY - visualisation size (default 800, 600)" - , " zoom - pixel replication factor (default 5)" - , " scale - feature size of visualisation (default 30)" - , " degree - number waves to sum for each point (default 5)" - , "" - , " You'll want to run this with +RTS -N to enable threads" ] - - -run :: Int -> Int -> Int -> Scale -> Degree -> IO () -run sizeX sizeY zoom scale degree - = animateField (InWindow "Crystal" (sizeX, sizeY) (10, 10)) - (zoom, zoom) - (quasicrystal scale degree) diff --git a/raster/Fluid/src-repa/Args.hs b/raster/Fluid/src-repa/Args.hs --- a/raster/Fluid/src-repa/Args.hs +++ b/raster/Fluid/src-repa/Args.hs @@ -1,230 +1,276 @@ -module Args where +module Args + ( parseArgs + , configDefault) +where import Config +import Model import Data.Array.Repa as R import Data.Array.Repa.Algorithms.Pixel as R import Data.Array.Repa.IO.BMP as R -import System.Console.GetOpt -import Data.IORef +import qualified Data.Vector.Unboxed as U import Prelude as P -import Control.Monad +import System.Exit +import Data.Char --- | Command line options. -loadConfig :: [String] -> IO Config -loadConfig args - = do - batchModeArg <- newIORef False - benchModeArg <- newIORef False - maxStepsArg <- newIORef 0 - widthArg <- newIORef 100 - itersArg <- newIORef 40 - scaleArg <- newIORef 5 - rateArg <- newIORef 25 - deltaArg <- newIORef 0.1 - diffArg <- newIORef 0.00001 - viscArg <- newIORef 0 - densArg <- newIORef 100 - velArg <- newIORef (20, 20) +parseArgs :: [String] -> Config -> IO Config +parseArgs args config + | [] <- args + = return config - densityBmpArg <- newIORef Nothing - velocityBmpArg <- newIORef Nothing - - let setWidthArg arg = writeIORef widthArg (read arg) - let setScaleArg arg = writeIORef scaleArg (read arg) - let setItersArg arg = writeIORef itersArg (read arg) - let setDeltaArg arg = writeIORef deltaArg (read arg) - let setDiffArg arg = writeIORef diffArg (read arg) - let setViscArg arg = writeIORef viscArg (read arg) - let setDensArg arg = writeIORef densArg (read arg) - let setVelArg arg = let a = read arg in writeIORef velArg (a, a) - let setRate arg = writeIORef rateArg (read arg) - let setMaxSteps arg = writeIORef maxStepsArg (read arg) - let setBatchArg = writeIORef batchModeArg True - let setBenchArg = writeIORef benchModeArg True - let setDensityBMP arg = writeIORef densityBmpArg (Just arg) - let setVelocityBMP arg = writeIORef velocityBmpArg (Just arg) + | "-batch" : rest <- args + = parseArgs rest + $ config { configBatchMode = True } - let options :: [OptDescr (IO ())] - options - = [Option [] ["batch"] (NoArg setBatchArg ) - "Run a fixed number of steps instead of displaying in a window.", + | "-frames" : path : rest <- args + = parseArgs rest + $ config { configFramesMode = Just path + , configBatchMode = True } - Option [] ["bench"] (NoArg setBenchArg ) - "Set standard initial conditions for benchmarking.", + | "-max" : num : rest <- args + , all isDigit num + = parseArgs rest + $ config { configMaxSteps = read num } - Option [] ["max"] (ReqArg setMaxSteps "INT") - "Quit after this number of steps.", + | "-size" : width : height : rest <- args + , all isDigit width + , all isDigit height + , width' <- read width + , height' <- read height + = parseArgs rest + $ config { configModelSize = (width', height') + , configInitialDensity = makeDensField_empty width' height' + , configInitialVelocity = makeVeloField_empty width' height' } - Option [] ["width"] (ReqArg setWidthArg "INT") - "Size of simulation (100)", + | "-unstable" : rest <- args + = parseArgs rest + $ config { configUnstable = True } - Option [] ["iters"] (ReqArg setItersArg "INT") - "Iterations for the linear solver (20)", + | "-iters" : num : rest <- args + , all isDigit num + = parseArgs rest + $ config { configIters = read num } - Option [] ["scale"] (ReqArg setScaleArg "INT") - "Width of cell in window (5)", + | "-scale" : int : rest <- args + = parseArgs rest + $ config { configScale = (read int, read int) } - Option [] ["rate"] (ReqArg setRate "INT") - "Frame rate of simulator (25)", + | "-rate" : int : rest <- args + = parseArgs rest + $ config { configRate = read int } - Option [] ["delta"] (ReqArg setDeltaArg "FLOAT") - "Length of time steps (0.1)", + | "-delta" : float : rest <- args + = parseArgs rest + $ config { configDelta = read float } - Option [] ["diffusion"] (ReqArg setDiffArg "FLOAT") - "Diffusion rate for the density (0)", + | "-diff" : float : rest <- args + = parseArgs rest + $ config { configDiff = read float } - Option [] ["viscosity"] (ReqArg setViscArg "FLOAT") - "Viscosity rate for the velocity (0)", + | "-diff-after" : int : rest <- args + , all isDigit int + = parseArgs rest + $ config { configDiffAfter = read int } - Option [] ["user-density"] (ReqArg setDensArg "FLOAT") - "Magnitude of a user inserted density (100)", + | "-visc" : float : rest <- args + = parseArgs rest + $ config { configVisc = read float } - Option [] ["user-velocity"] (ReqArg setVelArg "FLOAT") - "Magnitude of a user inserted velocity (20)", + | "-user-dens" : float : rest <- args + = parseArgs rest + $ config { configDensity = read float } - Option [] ["bmp-density"] (ReqArg setDensityBMP "FILE.bmp") - "File for initial fluid density", + | "-user-velo" : float : rest <- args + = parseArgs rest + $ config { configVisc = read float } - Option [] ["bmp-velocity"] (ReqArg setVelocityBMP "FILE.bmp") - "File for initial fluid velocity" - ] - case getOpt RequireOrder options args of - (actions, [], []) - -> foldl (>>) (return ()) actions + -- Initial Confditions ---------------------------------------------------- + | "-dens-bmp" : filePath : rest <- args + = do dens <- loadDensBMP filePath + let Z :. height :. width = extent dens + parseArgs rest + $ config { configInitialDensity = dens + , configModelSize = (width, height) } - ( _, nonOpts, []) - -> error $ "unrecognized arguments: " P.++ unwords nonOpts + | "-velo-bmp" : filePath : rest <- args + = do velo <- loadVeloBMP filePath + let Z :. height :. width = extent velo + parseArgs rest + $ config { configInitialVelocity = velo + , configModelSize = (width, height) } - ( _, _, msgs) - -> error $ concat msgs P.++ usageInfo "Usage: fluid [OPTION...]" options + | "-init-checks" : rest <- args + = do let (width, height) = configModelSize config + parseArgs rest + $ config { configInitialDensity = makeDensField_checks width height + , configInitialVelocity = makeVeloField_empty width height } + | "-init-man" : rest <- args + = do let (width, height) = configModelSize config + parseArgs rest + $ config { configInitialDensity = makeDensField_checks width height + , configInitialVelocity = makeVeloField_man width height } - batchMode <- readIORef batchModeArg - benchMode <- readIORef benchModeArg - maxSteps <- readIORef maxStepsArg - width <- readIORef widthArg - let height = width - iters <- readIORef itersArg - scale <- readIORef scaleArg - rate <- readIORef rateArg - delta <- readIORef deltaArg - diff <- readIORef diffArg - visc <- readIORef viscArg - dens <- readIORef densArg - vel <- readIORef velArg - densityBmp <- readIORef densityBmpArg - velocityBmp <- readIORef velocityBmpArg + | "-init-elk" : rest <- args + = do let (width, height) = configModelSize config + parseArgs rest + $ config { configInitialDensity = makeDensField_checks width height + , configInitialVelocity = makeVeloField_elk width height } + | otherwise + = do printUsage + exitWith ExitSuccess - -- Load the initial desity bmp if we were given one. - let mkInitialDensity - -- Load density from a .bmp, using the luminance as - -- the scalar density value. - | Just filePath <- densityBmp - = do result <- readImageFromBMP filePath - let arr = case result of - Right arr' -> arr' - Left err -> error $ show err - let Z :. height' :. width' - = extent arr +printUsage :: IO () +printUsage + = putStr + $ unlines + [ "gloss-fluid [flags]" + , " -batch Run a fixed number of steps instead of displaying in a window." + , " -frames Dump all frames to .bmp files (implies -batch)" + , " -max Quit after this number of steps." + , " -size Size of simulation. (100)" + , " -unstable Use the unstable linear solver (False)" + , " -iters Iterations for the linear solver. (40)" + , " -scale Width of a cell in the window. (5)" + , " -rate Frame rate. (30)" + , " -delta Length of time step. (0.1)" + , " -diff Diffusion rate for the density. (0)" + , " -diff-after Trigger diffusion after this step. (0)" + , " -visc Diffusion rate for the velocity. (0)" + , " -user-dens Magnitude of user inserted density. (100)" + , " -user-velo Magnitude of user inserted velocity. (20)" + , " -bmp-dens File for initial fluid density." + , " -bmp-velo File for initial fluid velocity." + , "" + , " Run this with +RTS -N -qa -qg to enable threads." + , "" ] - when (height /= height' || width /= width') - $ error "Fluid: bmp size does not match --width argument" - density <- computeUnboxedP - $ R.map floatLuminanceOfRGB8 arr - return density +configDefault :: Config +configDefault + = let modelW = 100 + modelH = 100 + in Config + { configRate = 30 + , configMaxSteps = 0 + , configBatchMode = False + , configFramesMode = Nothing + , configModelSize = (modelW, modelH) + , configScale = (5, 5) + , configUnstable = False + , configIters = 40 + , configDelta = 0.1 + , configDiff = 0 + , configDiffAfter = 0 + , configVisc = 0 + , configDensity = 100 + , configVelocity = (20, 20) + , configInitialDensity = makeDensField_empty modelW modelH + , configInitialVelocity = makeVeloField_empty modelW modelH } - | benchMode - = let width' = fromIntegral width - yc = fromIntegral (width `div` 2) - xc = fromIntegral (width `div` 2) - - in return - $ R.fromListUnboxed (Z :. height :. width) - $ [ let x' = fromIntegral (x - 1) - y' = fromIntegral (y - 1) - xk1 = cos (10 * (x' - xc) / width') - yk1 = cos (10 * (y' - yc) / width') - d1 = xk1 * yk1 - in if (d1 < 0) then 0 else d1 - | y <- [1..width] - , x <- [1..width] ] - -- No density file given, so just set the field to zero. - | otherwise - = return - $ R.fromListUnboxed (Z :. height :. width) - $ replicate (height * width) 0 +-- | Load a density field from a BMP file. +loadDensBMP :: FilePath -> IO DensityField +loadDensBMP filePath + = do result <- readImageFromBMP filePath + let arr = case result of + Right arr' -> arr' + Left err -> error $ show err + density <- computeUnboxedP + $ R.map floatLuminanceOfRGB8 arr - initialDensity <- mkInitialDensity + return density - -- Load the initial velocity bmp if we were given one - let mkInitialVelocity - -- Load - | Just filePath <- velocityBmp - = do result <- readImageFromBMP filePath - let arr = case result of - Right arr' -> arr' - Left err -> error $ show err +-- | Load velocity field from a BMP file. +loadVeloBMP :: FilePath -> IO VelocityField +loadVeloBMP filePath + = do result <- readImageFromBMP filePath + let arr = case result of + Right arr' -> arr' + Left err -> error $ show err - let Z :. height' :. width' - = extent arr + let {-# INLINE conv #-} + conv (r, g, _b) + = let r' = fromIntegral (-128 + fromIntegral r :: Int) + g' = fromIntegral (-128 + fromIntegral g :: Int) + in (r' * 0.0001, g' * 0.0001) - when (height /= height' || width /= width') - $ error "Fluid: bmp size does not match --width argument" + velocity <- computeUnboxedP $ R.map conv arr - let {-# INLINE conv #-} - conv (r, g, _b) - = let r' = fromIntegral (-128 + fromIntegral r :: Int) - g' = fromIntegral (-128 + fromIntegral g :: Int) - in (r' * 0.0001, g' * 0.0001) + return velocity - velocity <- computeUnboxedP $ R.map conv arr - return velocity - | benchMode - = let width' = fromIntegral width - yc = fromIntegral (width `div` 2) - xc = fromIntegral (width `div` 2) +------------------------------------------------------------------------------- +makeDensField_empty :: Int -> Int -> Array U DIM2 Float +makeDensField_empty width height + = R.fromUnboxed (Z :. height :. width) + $ U.replicate (width * height) 0 + + +makeDensField_checks :: Int -> Int -> DensityField +makeDensField_checks width height + = let height' = fromIntegral height + xc = fromIntegral (width `div` 2) + yc = fromIntegral (height `div` 2) - in return - $ R.fromListUnboxed (Z :. height :. width) - $ [ let x' = fromIntegral x - y' = fromIntegral y - xk2 = cos (15 * (x' - xc) / width') - yk2 = cos (15 * (y' - yc) / width') - d2 = xk2 * yk2 / 5 - in (0, d2) - | y <- [0..width-1] - , x <- [0..width-1] ] + in R.fromListUnboxed (Z :. height :. width) + $ [ let x' = fromIntegral (x - 1) + y' = fromIntegral (y - 1) + tx = 10 * (x' - xc) / height' + ty = 10 * (y' - yc) / height' + xk1 = if abs tx > 3*pi/2 then 0 else cos tx + yk1 = if abs ty > 3*pi/2 then 0 else cos ty + d1 = xk1 * yk1 + in if (d1 < 0) then 0 else d1 + | y <- [1 .. height] + , x <- [1 .. width] ] - -- No velocity file given, so just set the field to zero. - | otherwise - = return - $ R.fromListUnboxed (Z :. height :. width) - $ replicate (height * width) (0, 0) - initialVelocity <- mkInitialVelocity +------------------------------------------------------------------------------- +makeVeloField_empty :: Int -> Int -> Array U DIM2 (Float, Float) +makeVeloField_empty width height + = R.fromUnboxed (Z :. height :. width) + $ U.replicate (width * height) (0, 0) - return Config - { configRate = rate - , configWindowSize = (scale * width, scale * width) - , configMaxSteps = maxSteps - , configBatchMode = batchMode - , configModelSize = (width, width) - , configIters = iters - , configDelta = delta - , configDiffusion = diff - , configViscosity = visc - , configDensity = dens - , configVelocity = vel - , configInitialDensity = initialDensity - , configInitialVelocity = initialVelocity } +makeVeloField_man :: Int -> Int -> VelocityField +makeVeloField_man width height + = let height' = fromIntegral height + xc = fromIntegral (width `div` 2) + yc = fromIntegral (height `div` 2) + + in R.fromListUnboxed (Z :. height :. width) + $ [ let x' = fromIntegral x + y' = fromIntegral y + xk2 = cos (19 * (x' - xc) / height') + yk2 = cos (17 * (y' - yc) / height') + d2 = xk2 * yk2 / 5 + in (0, d2) + | y <- [0..height - 1] + , x <- [0..width - 1] ] + + +makeVeloField_elk :: Int -> Int -> VelocityField +makeVeloField_elk width height + = let height' = fromIntegral height + xc = fromIntegral (width `div` 2) + yc = fromIntegral (height `div` 2) + + in R.fromListUnboxed (Z :. height :. width) + $ [ let x' = fromIntegral x + y' = fromIntegral y + tx = 12 * (x' - xc) / height' + ty = 12 * (y' - yc) / height' + xk2 = cos tx + yk2 = -cos ty + d2 = xk2 * yk2 / 5 + in (0, d2) + | y <- [0 .. height - 1] + , x <- [0 .. width - 1] ] diff --git a/raster/Fluid/src-repa/Config.hs b/raster/Fluid/src-repa/Config.hs --- a/raster/Fluid/src-repa/Config.hs +++ b/raster/Fluid/src-repa/Config.hs @@ -1,18 +1,16 @@ module Config ( Config (..) - , configScale) + , configWindowSize) where import Model -import Data.Array.Repa +import Data.Array.Repa as R + data Config = Config - { -- | Size of window in pixels - configWindowSize :: (Int, Int) - - -- | Simulation rate (frames per second) - , configRate :: Int + { -- | Simulation rate (frames per second) + configRate :: Int -- | Maximum number of steps in simulation , configMaxSteps :: Int @@ -20,9 +18,18 @@ -- | Whether to run in batch-mode, non-interactively. , configBatchMode :: Bool + -- | Whether to dump all frames to .bmp files. + , configFramesMode :: Maybe FilePath + -- | Number of cells in model. , configModelSize :: (Int, Int) + -- | Window scale. + , configScale :: (Int, Int) + + -- | Use the unstable linear solver. + , configUnstable :: !Bool + -- | Number of iterations to use in the linear solver , configIters :: !Int @@ -30,10 +37,13 @@ , configDelta :: !Delta -- | Diffusion rate. - , configDiffusion :: !Float + , configDiff :: !Float + -- | Apply diffusion after this step number. + , configDiffAfter :: !Int + -- | Fluid viscosity. - , configViscosity :: !Float + , configVisc :: !Float -- | Magnitude of density to add with user interface. , configDensity :: !Float @@ -48,10 +58,9 @@ , configInitialVelocity :: !(Array U DIM2 (Float, Float)) } +configWindowSize :: Config -> (Int, Int) +configWindowSize config + = let (modelX, modelY) = configModelSize config + (scaleX, scaleY) = configScale config + in (modelX * scaleX, modelY * scaleY) -configScale :: Config -> (Float, Float) -configScale config - = let (windowWidth, windowHeight) = configWindowSize config - (modelWidth, modelHeight) = configModelSize config - in ( fromIntegral windowWidth / fromIntegral modelWidth - , fromIntegral windowHeight / fromIntegral modelHeight) diff --git a/raster/Fluid/src-repa/Main.hs b/raster/Fluid/src-repa/Main.hs --- a/raster/Fluid/src-repa/Main.hs +++ b/raster/Fluid/src-repa/Main.hs @@ -25,11 +25,11 @@ main :: IO () main = do -- Parse the command-line arguments. - args <- getArgs - config <- loadConfig args + args <- getArgs + config <- parseArgs args configDefault -- Setup the initial fluid model. - let model = initModel + let model = initModel (configInitialDensity config) (configInitialVelocity config) @@ -45,13 +45,14 @@ -- | Run the simulation interactively. runInteractive :: Config -> Model -> IO () runInteractive config model0 - = playIO (InWindow "Stam's stable fluid. Use left-click right-drag to add density / velocity." + = let (scaleX, scaleY) = configScale config + in playIO (InWindow "Stam's stable fluid. Use left-click right-drag to add density / velocity." (configWindowSize config) (20, 20)) black (configRate config) model0 - (pictureOfModel (configScale config)) + (pictureOfModel (fromIntegral scaleX, fromIntegral scaleY)) (\event model -> return $ userEvent config event model) (\_ -> stepFluid config) @@ -59,29 +60,31 @@ -- | Run in batch mode and dump a .bmp of the final state. runBatchMode :: Config -> Model -> IO () runBatchMode config model - | stepsPassed model >= configMaxSteps config - = return () + | stepsPassed model >= configMaxSteps config + = return () - | otherwise - = do model' <- stepFluid config model - runBatchMode config model' + | otherwise + = do model' <- stepFluid config model + case configFramesMode config of + Nothing -> return () + Just path + -> do putStrLn $ "frame " ++ show (stepsPassed model) + outputBMP path (stepsPassed model) (densityField model) + runBatchMode config model' + + -- Function to step simulator one step forward in time stepFluid :: Config -> Model -> IO Model -stepFluid config m@(Model df ds vf vs cl step cb) - | step >= configMaxSteps config - , configMaxSteps config > 0 - = return m +stepFluid config m@(Model df ds vf vs cl step dv cb) + | step >= configMaxSteps config + , configMaxSteps config > 0 + = return m - | otherwise - = do performGC + | otherwise + = do performGC vf' <- velocitySteps config step vf vs df' <- densitySteps config step df ds vf' - return $ Model df' Nothing vf' Nothing cl (step + 1) cb - - - - - + return $ Model df' Nothing vf' Nothing cl (step + 1) dv cb diff --git a/raster/Fluid/src-repa/Model.hs b/raster/Fluid/src-repa/Model.hs --- a/raster/Fluid/src-repa/Model.hs +++ b/raster/Fluid/src-repa/Model.hs @@ -70,6 +70,7 @@ , clickLoc :: Maybe (Int, Int) , stepsPassed :: Int , currButton :: CurrentButton + , drawVelocity :: Bool } @@ -91,7 +92,8 @@ , velocitySource = Nothing , clickLoc = Nothing , stepsPassed = 0 - , currButton = None } + , currButton = None + , drawVelocity = False } {-# INLINE initModel #-} @@ -99,7 +101,7 @@ -- | Function to convert the Model into a Bitmap for displaying in Gloss pictureOfModel :: Monad m => (Float, Float) -> Model -> m G.Picture pictureOfModel (scaleX, scaleY) m - = let (Z :. width' :. height') = R.extent $ densityField m + = let (Z :. height' :. width') = R.extent $ densityField m width = fromIntegral width' height = fromIntegral height' @@ -107,10 +109,32 @@ (arrDensity :: Array F DIM2 Word32) <- computeP $ R.map pixel32OfDensity $ densityField m - return $ G.Scale scaleX scaleY - $ G.bitmapOfForeignPtr width height + let picVel :: G.Picture + picVel + | drawVelocity m + = G.Translate (- width / 2) (- height / 2) + $ G.Color (G.light $ G.light G.red) + $ G.Pictures + [ G.Line [(xf, yf), (xf + vx', yf + vy')] + | x <- [0, 5 .. width' - 1] + , y <- [0, 5 .. height' - 1] + , let xf = fromIntegral x + , let yf = fromIntegral y + , let (vx0, vy0) = velocityField m R.! (Z :. y :. x) + , let vx' = vx0 * 100 + , let vy' = vy0 * 100 ] + + | otherwise + = G.blank + + let picDens :: G.Picture + picDens = G.bitmapOfForeignPtr width' height' (R.toForeignPtr $ unsafeCoerce arrDensity) False + + return $ G.Scale scaleX scaleY + $ G.pictures [picDens, picVel ] + {-# NOINLINE pictureOfModel #-} @@ -145,10 +169,11 @@ -- Dump ----------------------------------------------------------------------- -- Writes bitmap data to test batch-mode ran correctly -outputBMP :: DensityField -> IO () -outputBMP df - = do arr <- computeUnboxedP $ R.map pixel8OfDensity df - R.writeImageToBMP "./out.bmp" arr +outputBMP :: FilePath -> Int -> DensityField -> IO () +outputBMP path step df + = do arr <- computeUnboxedP $ R.map pixel8OfDensity df + let step' = replicate (6 - length (show step)) '0' P.++ show step + R.writeImageToBMP (path P.++ step' P.++ ".bmp") arr outputPPM :: Int -> String -> Float -> Array U DIM2 Float -> IO () diff --git a/raster/Fluid/src-repa/Solve/Density.hs b/raster/Fluid/src-repa/Solve/Density.hs --- a/raster/Fluid/src-repa/Solve/Density.hs +++ b/raster/Fluid/src-repa/Solve/Density.hs @@ -17,12 +17,22 @@ -> VelocityField -> IO DensityField -densitySteps config _step df ds vf +densitySteps config step df ds vf = {-# SCC "Solve.densitySteps" #-} do df1 <- addSources (configDelta config) (configDensity config) ds df - df2 <- diffusion (configIters config) (configDelta config) (configDiffusion config) + let diff = if configDiffAfter config /= 0 + && step >= configDiffAfter config + then 0.0005 + else configDiff config + + let diffSolver + = if configUnstable config + then DiffUnstable + else DiffStable (configIters config) + + df2 <- diffusion diffSolver (configDelta config) diff df1 df' <- advection (configDelta config) vf df2 diff --git a/raster/Fluid/src-repa/Solve/Velocity.hs b/raster/Fluid/src-repa/Solve/Velocity.hs --- a/raster/Fluid/src-repa/Solve/Velocity.hs +++ b/raster/Fluid/src-repa/Solve/Velocity.hs @@ -25,7 +25,8 @@ vf1 <- addSources (configDelta config) (configVelocity config) vs vf - vf2 <- diffusion (configIters config) (configDelta config) (configViscosity config) + let diffSolver = DiffStable (configIters config) + vf2 <- diffusion diffSolver (configDelta config) (configVisc config) vf1 -- vf3 <- setBoundary vf2 diff --git a/raster/Fluid/src-repa/Stage/Advection.hs b/raster/Fluid/src-repa/Stage/Advection.hs --- a/raster/Fluid/src-repa/Stage/Advection.hs +++ b/raster/Fluid/src-repa/Stage/Advection.hs @@ -48,8 +48,9 @@ (((d00 ~*~ t0) ~+~ (d01 ~*~ t1)) ~*~ s0) ~+~ (((d10 ~*~ t0) ~+~ (d11 ~*~ t1)) ~*~ s1) where - _ :. _ :. width' = R.extent velField - !width = fromIntegral width' + _ :. height' :. width' = R.extent velField + !width = fromIntegral width' + !height = fromIntegral height' -- helper values !dt0 = delta * width @@ -65,7 +66,7 @@ | otherwise = x !y' | y < -0.5 = -0.5 - | y > width + 0.5 = width + 0.5 + | y > height + 0.5 = height + 0.5 | otherwise = y -- calculate discrete locations surrounding point @@ -86,7 +87,7 @@ | ii < 0 = zero | ii >= width' = zero | jj < 0 = zero - | jj >= width' = zero + | jj >= height' = zero | otherwise = get ix -- grab values from grid surrounding advected point diff --git a/raster/Fluid/src-repa/Stage/Diffusion.hs b/raster/Fluid/src-repa/Stage/Diffusion.hs --- a/raster/Fluid/src-repa/Stage/Diffusion.hs +++ b/raster/Fluid/src-repa/Stage/Diffusion.hs @@ -1,6 +1,7 @@ {-# LANGUAGE BangPatterns #-} module Stage.Diffusion - (diffusion) + ( DiffSolver (..) + , diffusion) where import Model import FieldElt @@ -10,30 +11,40 @@ import Data.Vector.Unboxed +data DiffSolver + = DiffStable Int + | DiffUnstable + + -- | Diffuse a field at a certain rate. diffusion :: (FieldElt a, Num a, Elt a, Unbox a) - => Int + => DiffSolver -> Delta -> Rate -> Field a -> IO (Field a) -diffusion !iters !delta !rate field +diffusion !solver !delta !rate field = {-# SCC diffusion #-} field `deepSeqArray` let _ :. _ :. width' = R.extent field !width = fromIntegral width' - !a = delta * rate * width * width - !c = 1 + 4 * a - in do - linearSolver field field a c iters + in case solver of + DiffUnstable + -> let !a = delta * rate * width * width + in unstableSolver field field a + DiffStable iters + -> let !a = delta * rate * width * width + !c = 1 + 4 * a + in linearSolver field field a c iters + {-# SPECIALIZE diffusion - :: Int -> Delta -> Rate + :: DiffSolver -> Delta -> Rate -> Field Float -> IO (Field Float) #-} {-# SPECIALIZE diffusion - :: Int -> Delta -> Rate + :: DiffSolver -> Delta -> Rate -> Field (Float, Float) -> IO (Field (Float, Float)) #-} diff --git a/raster/Fluid/src-repa/Stage/Linear.hs b/raster/Fluid/src-repa/Stage/Linear.hs --- a/raster/Fluid/src-repa/Stage/Linear.hs +++ b/raster/Fluid/src-repa/Stage/Linear.hs @@ -1,6 +1,7 @@ {-# LANGUAGE FlexibleContexts, BangPatterns #-} module Stage.Linear - (linearSolver) + ( linearSolver + , unstableSolver) where import Model import FieldElt as E @@ -12,6 +13,7 @@ import Prelude as P +------------------------------------------------------------------------------- linearSolver :: (FieldElt a, Source U a, Unbox a, Elt a, Num a, Show a) => Field a -- ^ Original field. @@ -83,4 +85,64 @@ | i == 0, j == -1 = Just 1 | otherwise = Nothing {-# INLINE linearSolverCoeffs #-} + + +-- Unstable ------------------------------------------------------------------- +unstableSolver + :: (FieldElt a, Source U a, Unbox a, Elt a, Num a, Show a) + => Field a -- ^ Original field. + -> Field a -- ^ Current field. + -> Float + -> IO (Field a) + +unstableSolver !origField !curField !a + = origField `deepSeqArray` curField `deepSeqArray` + do + let {-# INLINE zipFunc #-} + zipFunc !orig !new + = orig ~+~ (new ~*~ a) + + newField + <- computeUnboxedP + $ R.szipWith zipFunc origField + $ mapStencil2 (BoundConst 0) unstableSolverStencil curField + + return newField + + +{-# SPECIALIZE unstableSolver + :: Field Float + -> Field Float + -> Float + -> IO (Field Float) #-} + +{-# SPECIALIZE unstableSolver + :: Field (Float, Float) + -> Field (Float, Float) + -> Float + -> IO (Field (Float, Float)) #-} + + +unstableSolverStencil + :: FieldElt a + => Stencil DIM2 a + +unstableSolverStencil + = StencilStatic (Z:.3:.3) E.zero + (\ix val acc -> + case unstableSolverCoeffs ix of + Nothing -> acc + Just coeff -> acc ~+~ (val ~*~ coeff)) +{-# INLINE unstableSolverStencil #-} + + +unstableSolverCoeffs :: DIM2 -> Maybe Float +unstableSolverCoeffs (Z:.j:.i) + | i == 1, j == 0 = Just 1 + | i == -1, j == 0 = Just 1 + | i == 0, j == 1 = Just 1 + | i == 0, j == -1 = Just 1 + | i == 0, j == 0 = Just (-4) + | otherwise = Nothing +{-# INLINE unstableSolverCoeffs #-} diff --git a/raster/Fluid/src-repa/Stage/Project.hs b/raster/Fluid/src-repa/Stage/Project.hs --- a/raster/Fluid/src-repa/Stage/Project.hs +++ b/raster/Fluid/src-repa/Stage/Project.hs @@ -14,18 +14,18 @@ project iters field = {-# SCC project #-} field `deepSeqArray` - do let _ :. _ :. width = extent field + do let _ :. height :. width = extent field divergence <- {-# SCC "project.genDiv" #-} computeUnboxedP - $ fromFunction (Z:. width :. width) (genDivergence width field) + $ fromFunction (Z:. height :. width) (genDivergence width height field) p <- {-# SCC "project.linearSolver" #-} linearSolver divergence divergence 1 4 iters f' <- {-# SCC "project.apply" #-} computeUnboxedP - $ unsafeTraverse field id (projectElem width p) + $ unsafeTraverse field id (projectElem width height p) return f' {-# NOINLINE project #-} @@ -35,32 +35,34 @@ -- create a mass-conserving field. projectElem :: Int -- ^ Width of model. + -> Int -- ^ Height of model. -> Field Float -> (DIM2 -> (Float, Float)) -- ^ Get data from the regular field. -> DIM2 -- ^ Compute the value at this point. -> (Float, Float) -projectElem !width !p !get !pos@(Z:.j:.i) - = get pos ~-~ (0.5 * width' * (p0 - p1), - 0.5 * width' * (p2 - p3)) +projectElem !width !height !p !get !pos@(Z:.j:.i) + = get pos ~-~ (0.5 * width' * (p0 - p1), + 0.5 * height' * (p2 - p3)) where - !width' = fromIntegral width - !p0 = useIf (i < width - 1) (p `unsafeIndex` (Z :. j :. i+1)) - !p1 = useIf (i > 0) (p `unsafeIndex` (Z :. j :. i-1)) - !p2 = useIf (j < width - 1) (p `unsafeIndex` (Z :. j+1 :. i )) - !p3 = useIf (j > 0) (p `unsafeIndex` (Z :. j-1 :. i )) + !width' = fromIntegral width + !height' = fromIntegral height + !p0 = useIf (i < width - 1) (p `unsafeIndex` (Z :. j :. i+1)) + !p1 = useIf (i > 0) (p `unsafeIndex` (Z :. j :. i-1)) + !p2 = useIf (j < height - 1) (p `unsafeIndex` (Z :. j+1 :. i )) + !p3 = useIf (j > 0) (p `unsafeIndex` (Z :. j-1 :. i )) {-# INLINE projectElem #-} -- | Get an approximation of the gradient at this point. -genDivergence :: Int -> VelocityField -> DIM2 -> Float -genDivergence !width !f (Z :. j :. i) +genDivergence :: Int -> Int -> VelocityField -> DIM2 -> Float +genDivergence !width !height !f (Z :. j :. i) = (-0.5 * ((u0 - u1) + (v0 - v1))) / fromIntegral width where - (u0, _) = useIf (i < width - 1) (f `unsafeIndex` (Z:. j :. i+1)) - (u1, _) = useIf (i > 0) (f `unsafeIndex` (Z:. j :. i-1)) - ( _, v0) = useIf (j < width - 1) (f `unsafeIndex` (Z:. j+1 :. i )) - ( _, v1) = useIf (j > 0) (f `unsafeIndex` (Z:. j-1 :. i )) + (u0, _) = useIf (i < width - 1) (f `unsafeIndex` (Z:. j :. i+1)) + (u1, _) = useIf (i > 0) (f `unsafeIndex` (Z:. j :. i-1)) + ( _, v0) = useIf (j < height - 1) (f `unsafeIndex` (Z:. j+1 :. i )) + ( _, v1) = useIf (j > 0) (f `unsafeIndex` (Z:. j-1 :. i )) {-# INLINE genDivergence #-} diff --git a/raster/Fluid/src-repa/UserEvent.hs b/raster/Fluid/src-repa/UserEvent.hs --- a/raster/Fluid/src-repa/UserEvent.hs +++ b/raster/Fluid/src-repa/UserEvent.hs @@ -11,109 +11,74 @@ userEvent :: Config -> Event -> Model -> Model userEvent config (EventKey key keyState mods (x, y)) - (Model df ds vf vs cl sp _cb) + model -- Add velocity --------------------------------------------- | MouseButton G.RightButton <- key , Down <- keyState , (x',y') <- windowToModel config (x,y) - = Model { densityField = df - , densitySource = ds - , velocityField = vf - , velocitySource = vs - , clickLoc = Just (x',y') - , stepsPassed = sp - , currButton = M.RightButton - } + = model { clickLoc = Just (x',y') + , currButton = M.RightButton } -- Accept shift-leftbutton for people with trackpads | MouseButton G.LeftButton <- key , Down <- keyState , Down <- shift mods , (x',y') <- windowToModel config (x,y) - = Model { densityField = df - , densitySource = ds - , velocityField = vf - , velocitySource = vs - , clickLoc = Just (x',y') - , stepsPassed = sp - , currButton = M.RightButton - } + = model { clickLoc = Just (x',y') + , currButton = M.RightButton } | MouseButton G.RightButton <- key , Up <- keyState - , Just (locX, locY) <- cl + , Just (locX, locY) <- clickLoc model , (x',y') <- windowToModel config (x,y) - = Model { densityField = df - , densitySource = ds - , velocityField = vf - , velocitySource = Just (SourceDensity (Z:.locY:.locX) + = model { velocitySource = Just (SourceDensity (Z:.locY:.locX) (fromIntegral (locX-x'),fromIntegral (locY-y'))) , clickLoc = Nothing - , stepsPassed = sp - , currButton = M.None - } - + , currButton = M.None } -- Add density ---------------------------------------------- | MouseButton G.LeftButton <- key , Down <- keyState , (x',y') <- windowToModel config (x, y) - = Model { densityField = df - , densitySource = Just (SourceDensity (Z:.y':.x') 1) - , velocityField = vf - , velocitySource = vs - , clickLoc = cl - , stepsPassed = sp - , currButton = M.LeftButton - } + = model { densitySource = Just (SourceDensity (Z:.y':.x') 1) + , currButton = M.LeftButton } | MouseButton G.LeftButton <- key , Up <- keyState - = Model { densityField = df - , densitySource = ds - , velocityField = vf - , velocitySource = vs - , clickLoc = cl - , stepsPassed = sp - , currButton = M.None - } - + = model { currButton = M.None } - -- Reset model ---------------------------------------------- + -- Reset model | Char 'r' <- key , Down <- keyState = initModel (configInitialDensity config) (configInitialVelocity config) - -- Quit program --------------------------------------------- + -- Toggle velocity display + | Char 'v' <- key + , Down <- keyState + = model { drawVelocity = not $ drawVelocity model } + + -- Quit program | Char 'q' <- key , Down <- keyState = error "Quitting" -userEvent config (EventMotion (x, y)) - (Model df _ds vf vs cl sp M.LeftButton) - | (x',y') <- windowToModel config (x, y) - = Model { densityField = df - , densitySource = Just (SourceDensity (Z:.y':.x') 1) - , velocityField = vf - , velocitySource = vs - , clickLoc = cl - , stepsPassed = sp +userEvent config (EventMotion (x, y)) model + | M.LeftButton <- currButton model + , (x',y') <- windowToModel config (x, y) + = model { densitySource = Just (SourceDensity (Z:.y':.x') 1) , currButton = M.LeftButton } -userEvent config (EventMotion (x, y)) - (Model df ds vf _vs (Just (clx, cly)) sp M.RightButton) - | (x', y') <- windowToModel config (x,y) - = Model { densityField = df - , densitySource = ds - , velocityField = vf - , velocitySource = Just (SourceDensity (Z:.y':.x') +userEvent config (EventMotion (x, y)) model + | (x', y') <- windowToModel config (x,y) + , Just (clx, cly) <- clickLoc model + , M.RightButton <- currButton model + = model { velocitySource = Just (SourceDensity (Z:.y':.x') (fromIntegral (clx-x'), fromIntegral (cly-y'))) , clickLoc = Just (x',y') - , stepsPassed = sp , currButton = M.RightButton } @@ -126,8 +91,8 @@ where (scaleX, scaleY) = configScale config (windowWidth, windowHeight) = configWindowSize config - x' = round ((x + (fromIntegral windowWidth / 2)) / scaleX) - y' = round ((y + (fromIntegral windowHeight / 2)) / scaleY) + x' = round ((x + (fromIntegral windowWidth / 2)) / fromIntegral scaleX) + y' = round ((y + (fromIntegral windowHeight / 2)) / fromIntegral scaleY)