diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c) 2012, Dmitry Dzhus + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Dmitry Dzhus nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/README b/README new file mode 100644 --- /dev/null +++ b/README @@ -0,0 +1,52 @@ +This package provides CLI to dsmc library to perform rarefied gas flow +simulation. The tool is named `dsmc-runner`. Another tool, +`dsmc-caster` is used to render CSG bodies to check if the body +geometry is defined correctly. + +Both programs can be run with `-?` flag to obtain usage informaiton +and help of command-line flags. + +To run programs using multiple threads, specify `+RTS -N` flag. + +Extra data files which are installed with this package can be found in +`share` directory used by Cabal. + +* dsmc-runner + +To perform a simulation, input parameters for DSMC must be defined in +«problem file», which uses INI-like format, see `example.prb` included +in the package. + +To define the body (which is then referenced from the problem file), +`.geo` format is used, see `body.geo`, `smiley.geo` and `reentry2.geo` +for example. + +The runner dumps macroscopic data in text format to the file specified +using `-m` option, which can be then analyzed using visualization +software like ParaView. + +To sum up, + + $ dsmc-runner example.prb -m macro.txt +RTS -N + + 2012-06-13 11:32:02 MSK: [DEBUG] Program started + 2012-06-13 11:32:02 MSK: [DEBUG] Using problem file: example.prb + 2012-06-13 11:32:02 MSK: [DEBUG] Macroscopic data output file: macro.txt + 2012-06-13 11:32:02 MSK: [DEBUG] Not writing ensemble data + 2012-06-13 11:32:02 MSK: [DEBUG] Using body definition file: ./reentry2.geo + 2012-06-13 11:32:02 MSK: [DEBUG] Starting simulation + 2012-06-13 11:32:04 MSK: [DEBUG] Not steady yet; particles count: 1662080 + 2012-06-13 11:32:05 MSK: [DEBUG] Not steady yet; particles count: 1661946 + 2012-06-13 11:32:06 MSK: [DEBUG] Not steady yet; particles count: 1662026 + … + +* dsmc-caster + +dsmc-caster is a program used to view bodies defined using `.geo` +format. + + $ dsmc-caster smiley.geo +RTS -N + +Drag with left mouse button to rotate the object and mouse wheel to +zoom in/out. More command-line options are available, see usage help +with `-?` flag. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/body.geo b/body.geo new file mode 100644 --- /dev/null +++ b/body.geo @@ -0,0 +1,8 @@ +solid cap = sphere (0, 0, -0.8; 3.1); +solid body = cone (0, 0, -5; 5; 0, 0, 0; 3); +solid bottom = sphere (0, 0, 0; 5); + +solid upper = cap or body; +solid whole = upper and bottom; + +tlo whole; diff --git a/dsmc-tools.cabal b/dsmc-tools.cabal new file mode 100644 --- /dev/null +++ b/dsmc-tools.cabal @@ -0,0 +1,49 @@ +name: dsmc-tools +description: CLI tools for dsmc library +version: 0.1.0.0 +synopsis: DSMC toolkit for rarefied gas dynamics +license: BSD3 +license-file: LICENSE +author: Dmitry Dzhus +maintainer: dima@dzhus.org +category: Physics + +data-files: README, example.prb, body.geo, smiley.geo, reentry2.geo + +build-type: Simple +cabal-version: >=1.8 +tested-with: GHC == 7.4.1 + +source-repository head + type: git + location: https://github.com/dzhus/dsmc-tools/ + +executable dsmc-caster + ghc-options: -Wall -O2 -funbox-strict-fields -Odph -rtsopts -threaded -fno-liberate-case -funfolding-use-threshold1000 -funfolding-keeness-factor1000 -fllvm -optlo-O3 + main-is: src/Caster.hs + + -- other-modules: + + build-depends: + base == 4.*, + cmdargs == 0.9.*, + dsmc == 0.1.*, + gloss == 1.7.*, + gloss-raster == 1.7.*, + strict == 0.3.* + +executable dsmc-runner + ghc-options: -Wall -O2 -funbox-strict-fields -Odph -rtsopts -threaded -fno-liberate-case -funfolding-use-threshold1000 -funfolding-keeness-factor1000 -fllvm -optlo-O3 + main-is: src/Runner.hs + + build-depends: + base == 4.*, + bytestring == 0.9.*, + cmdargs == 0.9.*, + ConfigFile == 1.1.*, + dsmc == 0.1.*, + hslogger == 1.2.*, + transformers == 0.3.*, + mtl == 2.1.*, + repa == 3.2.*, + vector == 0.9.* diff --git a/example.prb b/example.prb new file mode 100644 --- /dev/null +++ b/example.prb @@ -0,0 +1,52 @@ +[Body] +# Path to body definition file +definition = ./reentry2.geo +volume_points = 1000 + +[Flow] +# Absolute temperature, in K +t = 300 + +# Concentration +n = 1e20 + +# Molecular mass +m = 29 + +# Flow velocity components, in m/s +velocity = (1000, 0, 0) + +# How many real particles a single simulator particle represents +fn = 5e16 + +[Surface] +temperature = 300 + +[Macroscopic] +# Macroscopic grid sizes, in m +mx = 0.1 +my = 0.1 +mz = 0.1 + +[Domain] +# Domain dimensions by x, y, z in m +w = 12 +l = 12 +h = 6 + +[Simulation] +dt = 1e-5 + +# For how many steps in steady regime to collect macroscopic data +steady_steps = 500 + +# Consider the system steady when the ratio of total particle counts +# on two consecutive steps differs from unit for less than this number +steady_epsilon = 1e-6 + +# Interface domains extrusion +ex = 0.04 + +# Start with empty simulation domain if true, fill it with particles +# distributed according to flow parameters upon start otherwise +empty_start = false diff --git a/reentry2.geo b/reentry2.geo new file mode 100644 --- /dev/null +++ b/reentry2.geo @@ -0,0 +1,5 @@ +solid s = sphere (0, 0, 0; 2); +solid s2 = sphere (1, 0, 0; 2.8); +solid whole = s and s2; + +tlo whole; diff --git a/smiley.geo b/smiley.geo new file mode 100644 --- /dev/null +++ b/smiley.geo @@ -0,0 +1,14 @@ +solid head = sphere (0, 0, 0; 10); +solid cut = plane (1, 0, 0; 1, 0, 0) and plane (-1, 0, 0; -1, 0, 0); + + +solid eyeballs = sphere (0, -4.2, 3; 2) or sphere (0, 4.2, 3; 2); +solid eyecuts = not eyeballs; + + +solid smile = cylinder (-2, 0, 0; 2, 0, 0; 5) and not cylinder (-3, 0, 2; 3, 0, 2; 5); + +solid smilecut = head and cut and eyecuts; +solid whole = smilecut and not smile; + +tlo whole; diff --git a/src/Caster.hs b/src/Caster.hs new file mode 100644 --- /dev/null +++ b/src/Caster.hs @@ -0,0 +1,207 @@ +{-# LANGUAGE DeriveDataTypeable #-} +{-# LANGUAGE RecordWildCards #-} +{-# LANGUAGE BangPatterns #-} + +{-| + +Standalone raycasting runner for traceable objects. + +|-} + +import Prelude hiding (reverse) + +import GHC.Float + +import Graphics.Gloss.Data.Color +import Graphics.Gloss.Data.Display +import qualified Graphics.Gloss.Data.Point as G +import qualified Graphics.Gloss.Interface.Pure.Game as G +import Graphics.Gloss.Raster.Field hiding (Point) + +import System.Console.CmdArgs.Implicit + +import qualified Data.Strict.Maybe as S +import qualified Data.Strict.Tuple as S + +import DSMC.Particles +import DSMC.Traceables +import DSMC.Traceables.Parser +import DSMC.Util.Vector + + +data InteractionMode = None | Rotate | Pan + + +-- | World state with observation point parameters and event helpers. +-- Roll is always 0. +data World = World { dist :: Double + -- ^ Distance to origin. + , pitch :: Double + , yaw :: Double + -- ^ Yaw of camera as if it was in origin. + , target :: Point + -- ^ Where camera looks at. + , holdPoint :: Maybe (Float, Float) + -- ^ Point where mouse button was held down. + , mode :: InteractionMode + } + + +-- | Static command line options for caster. +data Options = Options + { bodyDef :: FilePath + , width :: Int + , height :: Int + , pixels :: Int + , brightRGBA :: (Float, Float, Float, Float) + -- ^ Color for bright surfaces parallel to view plane (RGBA) + , darkRGBA :: (Float, Float, Float, Float) + -- ^ Color for dark surfaces perpendicular to view plane (RGBA) + } + deriving (Data, Typeable) + + +type Ray = Particle + + +-- | Pixels in meter at unit distance. +scaleFactor :: Double +scaleFactor = 50.0 + + +-- | Initial world. +start :: World +start = World 5 0 0 origin Nothing None + + +uncurry4 :: (a -> b -> c -> d -> e) -> (a, b, c, d) -> e +uncurry4 f (a, b, c, d) = f a b c d + + +-- | Scale deltas between hold and release coordinates by this number. +dragFactor :: Double +dragFactor = pi / 1800 + + +-- | Change distance by this amount per one mouse wheel step. +zoomFactor :: Double +zoomFactor = 0.01 + + +-- | Handle mouse drag to change pitch & yaw and mouse wheel to zoom. +handleEvents :: G.Event -> World -> World +handleEvents e w = + case e of + G.EventKey (G.MouseButton G.LeftButton) G.Down _ c -> + w{holdPoint = Just c, mode = Rotate} + G.EventKey (G.MouseButton G.RightButton) G.Down _ c -> + w{holdPoint = Just c, mode = Pan} + G.EventKey (G.MouseButton _) G.Up _ _ -> + w{holdPoint = Nothing} + G.EventKey (G.MouseButton G.WheelDown) _ _ _ -> + w{dist = dist w + zoomFactor} + G.EventKey (G.MouseButton G.WheelUp) _ _ _ -> + w{dist = dist w - zoomFactor} + G.EventKey (G.Char 'r') G.Down _ _ -> + w{ target = origin + , yaw = 0 + , pitch = 0} + G.EventMotion p@(x, y) -> + case (holdPoint w) of + Nothing -> w + Just (u, v) -> + let + xdelta = (float2Double (x - u)) * dragFactor + ydelta = (float2Double (y - v)) * dragFactor + in + case (mode w) of + Rotate -> w{ holdPoint = Just p + , yaw = (yaw w) - xdelta + , pitch = (pitch w) + ydelta + } + Pan -> + let + !(_, sX, sY) = buildCartesian (yaw w) (pitch w) + in + w{ holdPoint = Just p + , target = (target w) <+> (sX .^ xdelta) <-> (sY .^ ydelta) + } + _ -> w + _ -> w + + +casterField :: Int + -- ^ Window width. + -> Int + -- ^ Window height. + -> Int + -- ^ Pixels per point. + -> Body + -- ^ Body to show. + -> Color + -- ^ Bright color. + -> Color + -- ^ Dark color. + -> IO () +casterField width height pixels body bright' dark' = + let + display = InWindow "dsmc-tools CSG raycaster" (width, height) (0, 0) + makePixel :: World -> G.Point -> Color + !wS = fromIntegral (width `div` 2) + !hS = fromIntegral (height `div` 2) + makePixel !w !(x, y) = + let + !d = dist w + !wScale = -(wS * d / scaleFactor) + !hScale = (hS * d / scaleFactor) + !(n, sX, sY) = buildCartesian (yaw w) (pitch w) + !p = n .^ (-d) <+> target w + ray :: Ray + !ray = ((p + <+> (sX .^ ((float2Double x) * wScale)) + <+> (sY .^ ((float2Double y) * hScale))), n) + + !hp = trace body ray + in + case hp of + (((S.:!:) (HitPoint _ (S.Just hn)) _):_) -> + mixColors factor (1 - factor) bright' dark' + where + factor = abs $ double2Float $ invert n .* hn + _ -> white + {-# INLINE makePixel #-} + in + playField display (pixels, pixels) 5 start makePixel + handleEvents + (flip const) + + +-- | Read body def and program arguments, run the actual caster on +-- success. +main :: IO () +main = + let + sample = Options + { bodyDef = def &= argPos 0 &= typ "BODY-FILE" + , width = 500 &= help "Window width" + , height = 500 &= help "Window height" + , pixels = 1 + &= help ("Number of pixels to draw per point," ++ + "in each dimension") + , brightRGBA = (0.9, 0.9, 0.9, 1) + &= help ("Color for surface parallel to view plane," ++ + " with each component within [0..1]") + &= typ "R,G,B,A" + , darkRGBA = (0, 0, 0, 1) + &= help "Color for surface perpendicular to view plane" + &= typ "R,G,B,A" + } + &= program "dsmc-caster" + in do + Options{..} <- cmdArgs $ sample + body <- parseBodyFile bodyDef + case body of + Right b -> casterField width height pixels b + (uncurry4 makeColor brightRGBA) + (uncurry4 makeColor darkRGBA) + Left e -> error $ "Problem when reading body definition: " ++ e diff --git a/src/Runner.hs b/src/Runner.hs new file mode 100644 --- /dev/null +++ b/src/Runner.hs @@ -0,0 +1,246 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE DeriveDataTypeable #-} +{-# LANGUAGE RecordWildCards #-} +{-# LANGUAGE ScopedTypeVariables #-} + +module Main + +where + +import Control.Monad.Error +import Control.Exception + +import Data.ByteString.Char8 as B + +import Data.ConfigFile + +import Data.Functor + +import Data.Vector.Unboxed as VU hiding ((++)) + +import System.Console.CmdArgs.Implicit + +import System.IO + +import System.Log +import System.Log.Logger +import System.Log.Formatter +import System.Log.Handler (setFormatter) +import System.Log.Handler.Simple + +import DSMC +import DSMC.Domain +import DSMC.Macroscopic +import DSMC.Particles +import DSMC.Surface +import DSMC.Traceables.Parser +import DSMC.Util.Constants +import DSMC.Util.Vector + +import Data.Array.Repa as R hiding ((++)) + +import GHC.Conc.Sync + + +simSection :: String +simSection = "Simulation" + + +flowSection :: String +flowSection = "Flow" + + +domainSection :: String +domainSection = "Domain" + + +macroSection :: String +macroSection = "Macroscopic" + +bodySection :: String +bodySection = "Body" + + +surfaceSection :: String +surfaceSection = "Surface" + + +-- | Command-line options for the runner. +data Options = Options + { probDef :: FilePath + , macroFile :: Maybe FilePath + , ensFile :: Maybe FilePath + } + deriving (Data, Typeable) + +main :: IO () +main = + let + sample = Options + { probDef = def &= argPos 0 &= typ "PROBLEM-FILE" + , macroFile = Nothing &= help "Macroscopic data output file" + , ensFile = Nothing &= help "Ensemble data output file" + } + &= program "dsmc-runner" + in do + -- | Setup logging to stderr + let lf = simpleLogFormatter "$time: [$prio] $msg" + lh <- (flip setFormatter) lf <$> verboseStreamHandler stderr DEBUG + updateGlobalLogger rootLoggerName (setLevel DEBUG . setHandlers [lh] ) + + debugM rootLoggerName "Program started" + + Options{..} <- cmdArgs $ sample + + debugM rootLoggerName ("Using problem file: " ++ probDef) + + debugM rootLoggerName $ + maybe "Not writing macroscopic data" + ("Macroscopic data output file: " ++) macroFile + + debugM rootLoggerName $ + maybe "Not writing ensemble data" + ("Ensemble data output file: " ++) ensFile + + -- Use number of cores as split factor for parallelizing + -- stochastic processes. + nsplit <- getNumCapabilities + + (hr :: Either IOException Handle) <- try (openFile probDef ReadMode) + case hr of + Left e -> errorM rootLoggerName $ show e + Right hdl -> do + res <- runErrorT $ do + cp <- join $ liftIO $ readhandle emptyCP hdl + -- Must refactor this somehow (actually, write CmdArgs + -- for .ini files) + n <- get cp flowSection "n" + t <- get cp flowSection "t" + m <- get cp flowSection "m" + v <- get cp flowSection "velocity" + sw <- get cp flowSection "fn" + w <- get cp domainSection "w" + l <- get cp domainSection "l" + h <- get cp domainSection "h" + ex <- get cp simSection "ex" + dt <- get cp simSection "dt" + ssteps <- get cp simSection "steady_steps" + sepsilon <- get cp simSection "steady_epsilon" + emptys <- get cp simSection "empty_start" + vpts <- get cp bodySection "volume_points" + mx <- get cp macroSection "mx" + my <- get cp macroSection "my" + mz <- get cp macroSection "mz" + bodyDef <- get cp bodySection "definition" + -- Interface to diffuse model only + surfTemp <- get cp surfaceSection "temperature" + return $ (Flow n t (m * amu) v sw, ex, dt, + ssteps, sepsilon, emptys, + vpts, + mx, my, mz, + makeDomain origin w l h, + bodyDef, + surfTemp) + case res of + Left e -> errorM rootLoggerName $ + "Problem when reading problem definition: " ++ show e + Right (flow, ex, dt, + ssteps, sepsilon, emptys, + vpts, + mx, my, mz, + domain, + bodyDef, + surfTemp) -> do + -- Try to read body from definition file + debugM rootLoggerName $ + "Using body definition file: " ++ bodyDef + body <- parseBodyFile bodyDef + case body of + Left e -> errorM rootLoggerName $ + "Problem when reading body definition: " ++ e + Right b -> + do + debugM rootLoggerName $ "Starting simulation" + -- Run the simulation, obtaining + -- iteration count and final + -- distributions + !(iters, e, macro) <- + simulate domain b flow dt + emptys ex sepsilon ssteps + (Diffuse surfTemp $ + DSMC.Particles.mass flow) + (mx, my, mz) vpts nsplit + debugM rootLoggerName "Simulation complete" + debugM rootLoggerName $ + "Particle count when finished: " ++ + (show $ ensembleSize e) + debugM rootLoggerName $ + "Total iterations done: " ++ + (show iters) + case macroFile of + Just fp -> do + debugM rootLoggerName $ + "Writing macroscopic data to file file: " ++ fp + saveMacroscopic macro fp + Nothing -> return () + case ensFile of + Just fp -> do + debugM rootLoggerName $ + "Writing ensemble data to file: " ++ fp + saveEnsemble e fp + Nothing -> return () + return () + + +fd :: Double -> ByteString +fd = pack . show +{-# INLINE fd #-} + + +ws :: ByteString +ws = B.singleton ' ' + + +nl :: ByteString +nl = B.singleton '\n' + + +-- | Save Repa array in text file, one element per line. +saveRepa :: (Unbox e) => + (e -> ByteString) + -- ^ Formatter function. + -> ByteString + -- ^ Header line + -> R.Array U sh e + -> FilePath + -> IO () +saveRepa f header array path = do + h <- openFile path WriteMode + hPut h $ B.concat [header, nl] + VU.forM_ (R.toUnboxed array) (hPut h . f) + + +formatParticle :: Particle -> ByteString +formatParticle ((x, y, z), (u, v, w)) = + B.concat [fd x, ws, fd y, ws, fd z, ws, fd u, ws, fd v, ws, fd w, nl] +{-# INLINE formatParticle #-} + + +formatMacroParameters :: (Point, IntensiveMacroParameters) -> ByteString +formatMacroParameters ((x, y, z), (n, (u, v, w), p, t)) = + if n == 0.0 then B.empty + else B.concat [fd x, ws, fd y, ws, fd z, ws, + fd n, ws, + fd u, ws, fd v, ws, fd w, ws, + fd p, ws, fd t, nl] +{-# INLINE formatMacroParameters #-} + + +saveMacroscopic :: MacroField -> FilePath -> IO () +saveMacroscopic = + saveRepa formatMacroParameters $ + B.pack "x y z number_density u_avg v_avg w_avg pressure trans_temp" + + +saveEnsemble :: Ensemble -> FilePath -> IO () +saveEnsemble = saveRepa formatParticle $ B.pack "x y z u v w"