dsmc (empty) → 0.1.0.0
raw patch · 17 files changed
+2541/−0 lines, 17 filesdep +attoparsecdep +basedep +bytestringsetup-changed
Dependencies added: attoparsec, base, bytestring, containers, entropy, hslogger, mwc-random, parallel, primitive, repa, strict, transformers, vector
Files
- LICENSE +30/−0
- Setup.hs +2/−0
- dsmc.cabal +65/−0
- src/Control/Parallel/Stochastic.hs +65/−0
- src/DSMC.hs +157/−0
- src/DSMC/Cells.hs +287/−0
- src/DSMC/Domain.hs +246/−0
- src/DSMC/Macroscopic.hs +254/−0
- src/DSMC/Motion.hs +75/−0
- src/DSMC/Particles.hs +98/−0
- src/DSMC/Surface.hs +103/−0
- src/DSMC/Traceables.hs +486/−0
- src/DSMC/Traceables/Parser.hs +288/−0
- src/DSMC/Util.hs +107/−0
- src/DSMC/Util/Constants.hs +31/−0
- src/DSMC/Util/Vector.hs +169/−0
- src/Data/Splittable.hs +78/−0
+ LICENSE view
@@ -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.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ dsmc.cabal view
@@ -0,0 +1,65 @@+name: dsmc+description: Direct Simulation Monte Carlo is the numerical+ used to model the behavior of rarefied gas flows.+ This implementation supports complex bodies+ defined using Constructive Solid Geometry, using+ uniform grids and ray-casting. Specular, diffuse+ and CLL gas-surface interaction models are+ provided. Macroscopic parameters of number+ density, absolute velocity, pressure and+ translational temperature are obtained as the+ result of the simulation. The library employs+ parallelism on all steps of the DSMC algorithm.+ See the dsmc-tools package for command-line+ interfaces to the library.+version: 0.1.0.0+synopsis: DSMC library for rarefied gas dynamics+license: BSD3+license-file: LICENSE+author: Dmitry Dzhus+maintainer: dima@dzhus.org+category: Physics+build-type: Simple+cabal-version: >=1.8+tested-with: GHC == 7.4.1++source-repository head+ type: git+ location: https://github.com/dzhus/dsmc/++library+ ghc-options: -Wall -O2 -funbox-strict-fields -Odph -rtsopts -fno-liberate-case -funfolding-use-threshold1000 -funfolding-keeness-factor1000 -fllvm -optlo-O3+ hs-source-dirs: src++ exposed-modules:+ DSMC,+ DSMC.Cells,+ DSMC.Domain,+ DSMC.Surface,+ DSMC.Traceables,+ DSMC.Traceables.Parser,+ DSMC.Macroscopic,+ DSMC.Motion,+ DSMC.Particles,+ DSMC.Util.Constants,+ DSMC.Util.Vector++ other-modules:+ Control.Parallel.Stochastic,+ Data.Splittable,+ DSMC.Util++ build-depends:+ attoparsec == 0.10.*,+ base == 4.*,+ bytestring == 0.9.*,+ containers == 0.4.*,+ entropy == 0.2.*,+ hslogger == 1.2.*,+ mwc-random >= 0.12.0.1,+ parallel == 3.2.*,+ primitive == 0.4.*,+ repa == 3.2.*,+ strict == 0.3.*,+ transformers == 0.3.*,+ vector == 0.9.*
+ src/Control/Parallel/Stochastic.hs view
@@ -0,0 +1,65 @@+{-# LANGUAGE Rank2Types #-}++{-|++Parallel stochastic sampling for 'mwc-random' package.++-}++module Control.Parallel.Stochastic+ ( purifyRandomST+ , ParallelSeeds+ , parMapST+ , splitParMapST+ )++where++import Control.Monad.ST+import Control.Parallel.Strategies++import System.Random.MWC++import Data.Splittable+++-- | Convert ST action with PRNG state into a pure function of seed.+purifyRandomST :: (forall s.GenST s -> ST s a) -> Seed -> (a, Seed)+purifyRandomST f seed = runST $ do+ g <- restore seed+ r <- f g+ g' <- save g+ return (r, g')+{-# INLINE purifyRandomST #-}+++type RandomFunction source result = (forall s.GenST s -> source -> ST s result)+++-- | 'parMap' with 'rpar' over list of data and initial seeds using ST+-- action which takes single PRNG state; produce list of results and+-- used seeds.+parMapST :: RandomFunction a b -> [(a, Seed)] -> [(b, Seed)]+parMapST f = parMap rpar (\(p, seed) -> purifyRandomST (`f` p) seed)+{-# INLINE parMapST #-}+++-- | Split the given source, process subsources in parallel, return+-- combined results and used seeds.+splitParMapST :: (Split source, Combine result) =>+ RandomFunction source result+ -> source+ -> ParallelSeeds+ -> (result, ParallelSeeds)+splitParMapST f wholeSource oldSeeds =+ let+ sources = (splitIn (length oldSeeds) wholeSource)+ (results, newSeeds) = unzip $ parMapST f $ zip sources oldSeeds+ in+ (combine results, newSeeds)+{-# INLINE splitParMapST #-}+++-- | List of seeds which preserve PRNG states between runs of parallel+-- stochastic process sampling.+type ParallelSeeds = [Seed]
+ src/DSMC.hs view
@@ -0,0 +1,157 @@+{-# LANGUAGE BangPatterns #-}++{-|++DSMC is an algorithm used for simulating rarefied gas flows.++You define the simulation domain, the body inside this domain, gas+flow parameters and several other options. DSMC iteratively models the+behaviour of gas molecules according to time and space decoupling+scheme for the Boltzmann equation. The result of simulation is a field+of macroscopic parameters across the simulation domain.++-}++module DSMC+ ( motion+ , simulate+ )++where++import Control.Monad+import Control.Monad.IO.Class+import Data.Functor++import System.Log.Logger++import Control.Parallel.Stochastic++import DSMC.Cells+import DSMC.Domain+import DSMC.Macroscopic+import DSMC.Motion+import DSMC.Particles+import DSMC.Surface hiding (mass)+import DSMC.Traceables hiding (trace)+import DSMC.Util+++-- | Perform DSMC simulation, return total iterations count, final+-- particle distribution and field of averaged macroscopic parameters.+--+-- This is an IO action since system entropy source is polled for+-- seeds.+simulate :: Domain+ -> Body+ -> Flow+ -> Time+ -- ^ Time step.+ -> Bool+ -- ^ If true, start with empty domain. Add initial particle+ -- distribution to the domain otherwise.+ -> Double+ -- ^ Source reservoir extrusion.+ -> Double+ -- ^ Steadiness epsilon.+ -> Int+ -- ^ Step count limit in steady regime.+ -> Surface+ -- ^ Model for surface of body.+ -> (Double, Double, Double)+ -- ^ Spatial steps in X, Y, Z of grid used for macroscopic+ -- parameter sampling.+ -> Int+ -- ^ Use that many test points to calculate volume of every+ -- cell wrt body. Depends on Knudsen number calculated from+ -- cell size.+ -> Int+ -- ^ Split Lagrangian step into that many independent+ -- parallel processes.+ -> IO (Int, Ensemble, MacroField)+simulate domain body flow+ dt emptyStart ex sepsilon ssteps+ surface+ (mx, my, mz) volumePoints gsplit =+ let+ -- Simulate evolution of the particle system for one time+ -- step, updating seeds used for sampling stochastic+ -- processes.+ evolve :: (Ensemble, ParallelSeeds, DomainSeeds)+ -> (Ensemble, ParallelSeeds, DomainSeeds)+ evolve (ens, gseeds, dseeds) =+ let+ -- Inject new particles+ (e, dseeds') = openBoundaryInjection dseeds domain ex flow ens++ -- Lagrangian step+ (e', gseeds') = motion gseeds body dt surface e++ -- Filter out particles which left the domain+ e'' = clipToDomain domain e'+ in+ (e'', gseeds', dseeds')++ macroSubdiv :: Grid+ macroSubdiv = UniformGrid domain mx my mz++ -- Check if two consecutive particle ensemble states+ -- correspond to steady regime.+ stabilized :: Ensemble -> Ensemble -> Bool+ stabilized ens prevEns =+ (abs $+ ((fromIntegral $ ensembleSize ens) /+ (fromIntegral $ ensembleSize prevEns) - 1)) < sepsilon++ -- Helper which actually runs simulation and collects+ -- macroscopic data until enough samples in steady state are+ -- collected.+ sim1 :: (Ensemble, ParallelSeeds, DomainSeeds)+ -> Bool+ -- ^ True if steady regime has been reached.+ -> Int+ -- ^ Iteration counter.+ -> MacroSamplingMonad (Int, Ensemble, MacroField)+ sim1 !oldState@(ens, _, _) steady n =+ let+ !newState@(ens', _, _) = evolve oldState+ !newSteady = steady || stabilized ens' ens+ in do+ !enough <- case steady of+ False -> return False+ True -> updateSamples ens'+ liftIO $ debugM rootLoggerName $+ (if steady+ then "Steady state"+ else "Not steady yet") +++ "; particles count: " +++ (show $ ensembleSize ens')+ case enough of+ False -> sim1 newState newSteady (n + 1)+ True -> do+ (Just field) <- getField (mass flow) (statWeight flow)+ return (n, ens', field)+ in do+ -- Global seeds+ gs <- replicateM gsplit $ randomSeed++ -- Interface domain seeds+ (s1:s2:s3:s4:s5:s6:_) <- replicateM 6 randomSeed++ -- Seeds for cell volume calculation+ vs <- replicateM gsplit $ randomSeed++ -- Possibly sample initial particle distribution+ startEnsemble <- if emptyStart+ then return emptyEnsemble+ else do+ -- Forget the initial sampling seed+ is <- randomSeed+ return $ fst $ initializeParticles domain flow body is++ -- Start the process+ fst <$> runMacroSampling+ (sim1+ (startEnsemble, gs, (s1, s2, s3, s4, s5, s6))+ False 0)+ vs macroSubdiv body volumePoints ssteps
+ src/DSMC/Cells.hs view
@@ -0,0 +1,287 @@+{-# LANGUAGE BangPatterns #-}++{-|++Particle tracking for spatial grid for DSMC.++This module is used to sort (classify) particles into ordered vector+of cells for collision step or macroscopic parameter sampling. We do+not provide any special cell datatype since it varies which cell data+is required on every step, so only particles in every cell are stored.++Monad is provided for storing grid options during the whole program+run.++-}++module DSMC.Cells+ ( -- * Generic functions for cells+ Cells+ , CellContents+ , getCell+ , cellMap+ -- * Particle tracking+ , Classifier+ , classifyParticles+ -- * Grids+ , Grid(..)+ , makeUniformClassifier+ , makeUniformIndexer+ -- * Monadic interface+ , GridMonad+ , GridWares(..)+ , runGrid+ , cellVolumes+ )++where++import Prelude hiding (Just, Nothing, Maybe)++import Control.Monad.ST+import Control.Monad.Trans.Reader++import Data.Strict.Maybe++import qualified Data.Array.Repa as R+import qualified Data.Array.Repa.Repr.Vector as R++import qualified Data.Vector.Unboxed as VU+import qualified Data.Vector.Unboxed.Mutable as VUM+import qualified Data.Vector as V++import Control.Parallel.Stochastic++import DSMC.Domain+import DSMC.Particles+import DSMC.Traceables+import DSMC.Util+import DSMC.Util.Vector+++-- | Cell contents with particles.+type CellContents = VU.Vector Particle+++-- | Particles sorted by cells.+--+-- We store contents of all cells in a single densely packed unboxed+-- vector. Additionally cell count, cell starting positions in vector+-- (@s@) and cell sizes (@l@) are stored.+--+-- > s1 s2 s3+-- > | | |+-- > {[ooooooooo][oooo][oooooo]...}+-- > cell1 c2 c3+-- > l1=9 l2=4 l3=6+--+-- By using this storage scheme we allow fast particle classification+-- (see 'classifyParticles'). Slicing from contiguous memory block to+-- obtain contents of single cell is O(1) operation, and we maintain+-- data locality.+--+-- Note that any extra data about cells (like position or volume)+-- should be maintained separately from cell contents. We use this+-- approach because collision sampling and macroscopic parameter+-- calculation require different extra arguments.+data Cells = Cells !CellContents !Int !(VU.Vector Int) !(VU.Vector Int)+++-- | Assuming there's a linear ordering on all cells, Classifier must+-- yield index of cell for given particle.+type Classifier = Particle -> Int+++-- | Fetch contents of n-th cell.+getCell :: Cells+ -> Int+ -- ^ Cell index.+ -> Maybe CellContents+getCell !(Cells ens _ starts lengths) !n =+ case (lengths VU.! n) of+ 0 -> Nothing+ cl -> Just $ VU.slice (starts VU.! n) cl ens+{-# INLINE getCell #-}+++-- | Map a function over cell indices and contents of every cell.+cellMap :: (Int -> Maybe CellContents -> a) -> Cells -> R.Array R.D R.DIM1 a+cellMap f !cells@(Cells _ l _ _) =+ R.fromFunction (R.ix1 $ l) (\(R.Z R.:. cellNumber) ->+ f cellNumber $! getCell cells cellNumber)+++-- | Calculate cell numbers for particle ensemble.+classifyAll :: Classifier -> Ensemble -> (VU.Vector Int)+classifyAll classify ens = runST $ do+ classes' <- R.computeP $ R.map classify ens+ return $! R.toUnboxed classes'+++-- | Classify particle ensemble into @N@ cells using the classifier+-- function.+--+-- Classifier's extent must match @N@, yielding numbers between @0@+-- and @N-1@.+classifyParticles :: (Int, Classifier)+ -- ^ Cell count and classifier.+ -> Ensemble+ -> Cells+classifyParticles (cellCount, classify) ens' = runST $ do+ let ens = R.toUnboxed ens'+ particleCount = VU.length ens+ classes = classifyAll classify ens'++ -- Sequentially calculate particle indices inside cells and cell+ -- sizes.+ posns' <- VUM.replicate particleCount 0+ lengths' <- VUM.replicate cellCount 0+ iforM_ classes (\(particleNumber, cellNumber) -> do+ -- Increment cell particle count+ pos <- VUM.unsafeRead lengths' cellNumber+ VUM.unsafeWrite posns' particleNumber pos+ VUM.unsafeWrite lengths' cellNumber (pos + 1)+ return ())++ posns <- VU.unsafeFreeze posns'+ lengths <- VU.unsafeFreeze lengths'++ -- Starting positions for cells inside cell array+ let !starts = VU.prescanl' (+) 0 lengths++ -- Calculate indices for particles inside classified grand vector of+ -- cell contents (inverse mapping index)+ classifiedIds' <- VUM.replicate particleCount 0+ iforM_ classes (\(particleNumber, cellNumber) -> do+ let i = (starts VU.! cellNumber) + (posns VU.! particleNumber)+ VUM.unsafeWrite classifiedIds' i particleNumber+ return ())+ classifiedIds <- VU.unsafeFreeze classifiedIds'++ -- Fill the resulting array in parallel+ classifiedEns <- R.computeP $+ R.fromFunction+ (R.ix1 $ particleCount)+ (\(R.Z R.:. position) ->+ ens VU.! (classifiedIds VU.! position))++ return $! Cells (R.toUnboxed classifiedEns) cellCount starts lengths+++-- | Domain divided in uniform grid with given steps by X, Y and Z+-- axes.+data Grid = UniformGrid !Domain !Double !Double !Double+ deriving Show+++-- | Return grid cell count and classifier for a grid.+makeUniformClassifier :: Grid -> (Int, Classifier)+makeUniformClassifier (UniformGrid d@(Domain xmin _ ymin _ zmin _) hx hy hz) =+ (xsteps * ysteps * zsteps, classify)+ where+ (w, l, h) = getDimensions d+ xsteps = ceiling $ w / hx+ ysteps = ceiling $ l / hy+ zsteps = ceiling $ h / hz+ classify ((x, y, z), _) =+ let+ nx = floor $ (x - xmin) / hx+ ny = floor $ (y - ymin) / hy+ nz = floor $ (z - zmin) / hz+ in+ nx + ny * xsteps + nz * xsteps * ysteps+++-- | Function which maps cell numbers to central points of uniform+-- cells.+type Indexer = Int -> Point+++-- | Return indexer for a grid.+makeUniformIndexer :: Grid -> Indexer+makeUniformIndexer (UniformGrid d@(Domain xmin _ ymin _ zmin _) hx hy hz) =+ indefy+ where+ (w, l, _) = getDimensions d+ xsteps = ceiling $ w / hx+ ysteps = ceiling $ l / hy+ zf = xsteps * ysteps++ indefy i =+ let+ (nz, i') = i `divMod` zf+ z = zmin + fromIntegral nz * hz + hz / 2++ (ny, nx) = i' `divMod` ysteps+ y = ymin + fromIntegral ny * hy + hy / 2+ x = xmin + fromIntegral nx * hx + hx / 2+ in+ (x, y, z)+++-- | Build vector of domains corresponding to cells of grid.+gridDomains :: Grid -> V.Vector Domain+gridDomains g@(UniformGrid _ hx hy hz) =+ let+ ixer = makeUniformIndexer g+ (count, _) = makeUniformClassifier g+ in runST $ do+ doms <- R.computeP $ R.fromFunction (R.ix1 $ count)+ (\(R.Z R.:. cellNumber) ->+ makeDomain (ixer cellNumber) hx hy hz)+ return $ R.toVector doms+++-- | Calculate volumes of grid cells wrt body within the domain. For+-- every cell, 'freeVolume' is called with the domain of cell.+-- Calculation is performed in parallel.+--+-- Since our grid are static, this is usually done only once when the+-- grid is first defined. We throw away the used seeds.+cellVolumes :: ParallelSeeds+ -- ^ One-use seeds for cut cell volume approximation.+ -> Grid + -> Body + -> Int+ -> (VU.Vector Double)+cellVolumes seeds grid b testPoints =+ fst $+ splitParMapST (freeVolumes b testPoints)+ (gridDomains grid) seeds+++-- | Monad used to keep grid options and cell volumes. Due to the+-- low-level 'Cells' structure we use to store particles sorted in+-- cells, things may break badly if improper/inconsistent+-- classifier/indexer parameters are used with cells structure. It+-- also helps to maintain precalculated cell volumes. See+-- 'MacroSamplingMonad'.+type GridMonad = ReaderT GridWares DSMCRootMonad+++-- | Data stored in 'GridMonad'.+data GridWares =+ GridWares { classifier :: (Int, Classifier)+ -- ^ Cell count and classifier function.+ , indexer :: Int -> Point+ , volumes :: !(VU.Vector Double)+ -- ^ Vector of cell volumes.+ }+++-- | Run action using spatial subdivision.+runGrid :: GridMonad a + -> ParallelSeeds + -- ^ One-use seeds used for+ -> Grid+ -> Body+ -- ^ Body within the domain of the grid.+ -> Int+ -- ^ Use that many points to approximate every cell volume.+ -> DSMCRootMonad a+runGrid r seeds grid b testPoints =+ runReaderT r $+ GridWares+ (makeUniformClassifier grid)+ (makeUniformIndexer grid)+ (cellVolumes seeds grid b testPoints)
+ src/DSMC/Domain.hs view
@@ -0,0 +1,246 @@+{-# LANGUAGE BangPatterns #-}++{-|++Domain operations: defining domains; free flow boundary conditions &+clipping for DSMC steps.++PRNG required to sample molecular velocities implies monadic interface+for most of operations. We use functions specifically typed for 'ST'.++-}++module DSMC.Domain+ ( Domain(..)+ , getDimensions+ , getCenter+ , makeDomain+ -- * Flow boundary+ , initializeParticles+ , openBoundaryInjection+ , DomainSeeds+ , clipToDomain+ -- * Free volume calculation+ , freeVolume+ , freeVolumes+ )++where++import Control.Monad.ST++import qualified Data.Array.Repa as R+import qualified Data.Vector.Unboxed as VU+import qualified Data.Vector as V++import System.Random.MWC+import System.Random.MWC.Distributions (normal)++import Control.Parallel.Stochastic++import DSMC.Particles+import DSMC.Traceables+import DSMC.Util+import DSMC.Util.Constants+import DSMC.Util.Vector+++-- | Domain in which particles are spawned or system evolution is+-- simulated.+data Domain = Domain !Double !Double !Double !Double !Double !Double+ -- ^ Rectangular volume, given by min/max value on every+ -- dimension.+ deriving Show+++-- | Create a rectangular domain with center in the given point and+-- dimensions.+makeDomain :: Point+ -- ^ Center point.+ -> Double+ -- ^ X dimension.+ -> Double+ -- ^ Y dimension.+ -> Double+ -- ^ Z dimension.+ -> Domain+makeDomain !(x, y, z) !w !l !h =+ let+ xmin = x - w / 2+ ymin = y - l / 2+ zmin = z - h / 2+ xmax = x + w / 2+ ymax = y + l / 2+ zmax = z + h / 2+ in+ Domain xmin xmax ymin ymax zmin zmax+{-# INLINE makeDomain #-}+++-- | PRNG seeds used by particle generators.+type DomainSeeds = (Seed, Seed, Seed, Seed, Seed, Seed)+++-- | Calculate width, length and height of a domain, which are+-- dimensions measured by x, y and z axes, respectively.+getDimensions :: Domain -> (Double, Double, Double)+getDimensions (Domain xmin xmax ymin ymax zmin zmax) =+ (xmax - xmin, ymax - ymin, zmax - zmin)+{-# INLINE getDimensions #-}+++-- | Calculate geometric center of a domain.+getCenter :: Domain -> Point+getCenter (Domain xmin xmax ymin ymax zmin zmax) =+ (xmin + (xmax - xmin) / 2, ymin + (ymax - ymin) / 2, zmin + (zmax - zmin) / 2)+{-# INLINE getCenter #-}+++-- | Volume of domain.+volume :: Domain -> Double+volume !(Domain xmin xmax ymin ymax zmin zmax) =+ (xmax - xmin) * (ymax - ymin) * (zmax - zmin)+{-# INLINE volume #-}+++-- | Sample new particles inside a domain.+--+-- PRNG state implies this to be a monadic action.+spawnParticles :: Domain+ -> Flow+ -> GenST s+ -> ST s (VU.Vector Particle)+spawnParticles d@(Domain xmin xmax ymin ymax zmin zmax) flow g =+ let+ s = sqrt $ boltzmann * (temperature flow) / (mass flow)+ (u0, v0, w0) = velocity flow+ count = round $ (modelConcentration flow) * (volume d)+ in do+ VU.replicateM count $ do+ u <- normal u0 s g+ v <- normal v0 s g+ w <- normal w0 s g+ x <- uniformR (xmin, xmax) g+ y <- uniformR (ymin, ymax) g+ z <- uniformR (zmin, zmax) g+ return $ ((x, y, z), (u, v, w))+++-- | Pure version of 'spawnParticles'.+pureSpawnParticles :: Domain+ -> Flow+ -> Seed+ -> (VU.Vector Particle, Seed)+pureSpawnParticles d flow s = purifyRandomST (spawnParticles d flow) s+++-- | Fill the domain with particles for given flow parameters.+-- Particles inside the body are removed.+initializeParticles :: Domain+ -> Flow+ -> Body+ -> Seed+ -> (Ensemble, Seed)+initializeParticles d flow body s =+ let+ (res, s') = pureSpawnParticles d flow s+ ens = fromUnboxed1 res+ in + (filterEnsemble (not . inside body) ens, s')+++-- | Sample new particles in 6 interface domains along each side of+-- rectangular simulation domain and add them to existing ensemble.+--+-- This function implements open boundary condition for+-- three-dimensional simulation domain.+--+-- Interface domains are built on faces of simulation domain using+-- extrusion along the outward normal of the face.+--+-- In 2D projection:+--+-- > +-----------------++-- > | Interface1 |+-- > +--+-----------------+--++-- > |I3| Simulation |I4|+-- > | | domain | |+-- > +--+-----------------+--++-- > | I2 |+-- > +-----------------++--+-- Particles in every interface domain are spawned in parallel using+-- Strategies.+openBoundaryInjection :: DomainSeeds+ -> Domain+ -- ^ Simulation domain.+ -> Double+ -- ^ Interface domain extrusion length.+ -> Flow+ -> Ensemble+ -> (Ensemble, DomainSeeds)+openBoundaryInjection (s1, s2, s3, s4, s5, s6) domain ex flow ens =+ let+ (w, l, h) = getDimensions domain+ (cx, cy, cz) = getCenter domain+ d1 = makeDomain (cx - (w + ex) / 2, cy, cz) ex l h+ d2 = makeDomain (cx + (w + ex) / 2, cy, cz) ex l h+ d3 = makeDomain (cx, cy + (l + ex) / 2, cz) w ex h+ d4 = makeDomain (cx, cy - (l + ex) / 2, cz) w ex h+ d5 = makeDomain (cx, cy, cz - (h + ex) / 2) w l ex+ d6 = makeDomain (cx, cy, cz + (h + ex) / 2) w l ex+ v = [R.toUnboxed ens]+ (new, (s1':s2':s3':s4':s5':s6':_)) =+ unzip $+ parMapST (\g d -> spawnParticles d flow g) $+ zip [d1, d2, d3, d4, d5, d6] [s1, s2, s3, s4, s5, s6]+ in+ (fromUnboxed1 $ VU.concat (new ++ v), (s1', s2', s3', s4', s5', s6'))+++-- | Filter out particles which are outside of the domain.+clipToDomain :: Domain -> Ensemble -> Ensemble+clipToDomain (Domain xmin xmax ymin ymax zmin zmax) ens =+ let+ -- | Check if particle is in the domain.+ pred' :: Particle -> Bool+ pred' !((x, y, z), _) =+ xmax >= x && x >= xmin &&+ ymax >= y && y >= ymin &&+ zmax >= z && z >= zmin+ {-# INLINE pred' #-}+ in + filterEnsemble pred' ens+++-- | Volume of a domain unoccupied by a given body, in m^3.+--+-- We use Monte Carlo method to calculate the approximate body volume+-- and then subtract it from the overall domain volume.+freeVolume :: Domain + -> Body + -> Int + -- ^ Use that many points to approximate the body volume.+ -> GenST s+ -> ST s (Double)+freeVolume d@(Domain xmin xmax ymin ymax zmin zmax) body testPoints g = do+ points <- VU.replicateM testPoints $ do+ x <- uniformR (xmin, xmax) g+ y <- uniformR (ymin, ymax) g+ z <- uniformR (zmin, zmax) g+ return $ inside body ((x, y, z), (0, 0, 0))+ let occupiedPoints = VU.length $ VU.filter id points+ return $ (volume d) * + (fromIntegral (testPoints - occupiedPoints)) /+ (fromIntegral testPoints)+++-- | Sequential 'freeVolume' for a vector of domains.+freeVolumes :: Body+ -> Int+ -> GenST s+ -> V.Vector Domain+ -> ST s (VU.Vector Double)+freeVolumes body testPoints g doms =+ VU.generateM (V.length doms)+ (\i -> freeVolume (doms V.! i) body testPoints g)
+ src/DSMC/Macroscopic.hs view
@@ -0,0 +1,254 @@+{-# LANGUAGE BangPatterns #-}++{-|++Macroscopic parameters calculation.++We use regular spatial grid and time averaging for sampling. Sampling+should start after particle system has reached steady state. Samples+are then collected in each cell for a certain number of time steps.++Sampling is performed in 'MacroSamplingMonad' to ensure consistency of+averaging process. During sampling, basic parameters are calculated+like number of molecules per cell or mean square of thermal velocity.+After sampling these are used to derive final (intensive) parameters+like number density or temperature.++-}++module DSMC.Macroscopic+ ( MacroSamples+ , MacroField+ , BasicMacroParameters+ , IntensiveMacroParameters+ -- * Macroscopic sampling monad+ , MacroSamplingMonad+ , SamplingState(..)+ , runMacroSampling+ , updateSamples+ , getField+ )++where++import Control.Monad.Trans.Class (lift)+import Control.Monad.Trans.Reader+import Control.Monad.Trans.State.Strict++import qualified Data.Strict.Maybe as S++import qualified Data.Array.Repa as R+import qualified Data.Vector.Unboxed as VU++import Control.Parallel.Stochastic++import DSMC.Cells+import DSMC.Particles+import DSMC.Traceables+import DSMC.Util+import DSMC.Util.Constants+import DSMC.Util.Vector+++-- | Basic macroscopic parameters calculated in every cell: particle+-- count, mean absolute velocity, mean square of thermal velocity.+--+-- Particle count is non-integer because of averaging.+--+-- These are then post-processed into number density, flow velocity,+-- pressure and translational temperature.+--+-- Note the lack of root on thermal velocity!+type BasicMacroParameters = (Double, Vec3, Double)+++-- | Intensive macroscopic parameters available after averaging has+-- completed. These are: number density, absolute velocity, pressure+-- and translational temperature.+type IntensiveMacroParameters = (Double, Vec3, Double, Double)+++-- | Vector which stores averaged macroscropic parameters in each+-- cell.+--+-- If samples are collected for M iterations, then this vector is+-- built as a sum of vectors @V1, .. VM@, where @Vi@ is vector of+-- parameters sampled on @i@-th time step divided by @M@.+type MacroSamples = R.Array R.U R.DIM1 BasicMacroParameters+++-- | Array of central points of grid cells with averaged macroscopic+-- parameters attached to every point.+type MacroField = R.Array R.U R.DIM1 (Point, IntensiveMacroParameters)+++-- | Monad which keeps track of sampling process data and stores+-- options of macroscopic sampling.+--+-- GridMonad is used to ensure that only safe values for cell count+-- and classifier are used in 'updateSamples' and 'averageSamples'+-- (that may otherwise cause unbounded access errors). Note that+-- steady condition is not handled by this monad (instead, caller code+-- should decide when to start averaging).+--+-- Inner Reader Monad stores averaging steps setting.+type MacroSamplingMonad =+ StateT SamplingState (ReaderT Int GridMonad)+++-- | State of sampling process.+data SamplingState = None+ -- ^ Sampling has not started yet.+ | Incomplete Int MacroSamples+ -- ^ Sampling is in progress, not enough samples+ -- yet. Integer field indicates how many steps are+ -- left.+ | Complete MacroSamples+ -- ^ Averaging is complete, use 'getField' to+ -- unload the samples.+++makeIntensive :: Double+ -- ^ Mass of molecule.+ -> Double+ -- ^ Statistical weight of a simulator particle.+ -> Double+ -- ^ Cell volume.+ -> BasicMacroParameters+ -> IntensiveMacroParameters+makeIntensive !m !w !vol !(n, vel, c) =+ if (n == 0 || vol == 0)+ then (0, (0, 0, 0), 0, 0)+ else (numDens, vel, c * dens / 3, m * c / (3 * boltzmann))+ where+ numDens = n / vol * w+ dens = numDens * m+++-- | Fetch macroscopic field of intensive parameters if averaging is+-- complete.+getField :: Double+ -- ^ Mass of molecule.+ -> Double+ -- ^ Statistical weight of single molecule.+ -> MacroSamplingMonad (Maybe MacroField)+getField m w = do+ (cellCount, _) <- lift $ lift $ asks classifier+ ixer <- lift $ lift $ asks indexer+ vols <- lift $ lift $ asks volumes+ res <- get+ case res of+ Complete samples -> do+ let centralPoints = R.fromFunction (R.ix1 $ cellCount)+ (\(R.Z R.:. cellNumber) -> ixer cellNumber)+ realSamples = R.zipWith+ (makeIntensive m w)+ (fromUnboxed1 vols)+ samples+ f <- R.computeP $ R.zipWith (,) centralPoints realSamples+ return $ Just f+ _ -> return $ Nothing+++-- | Parameters in empty cell.+emptySample :: BasicMacroParameters+emptySample = (0, (0, 0, 0), 0)+++-- | Run 'MacroSamplingMonad' action with given sampling options and+-- return final 'Complete' state with macroscopic samples.+runMacroSampling :: MacroSamplingMonad r+ -> ParallelSeeds+ -> Grid+ -- ^ Grid used to sample macroscopic parameters.+ -> Body+ -> Int+ -- ^ Use that many points to approximate every cell volume.+ -> Int+ -- ^ Averaging steps count.+ -> DSMCRootMonad (r, SamplingState)+runMacroSampling f seeds grid body testPoints ssteps = + runGrid (runReaderT (runStateT f None) ssteps) seeds grid body testPoints+++-- | Create empty 'MacroSamples' array.+initializeSamples :: Int+ -- ^ Cell count.+ -> MacroSamples+initializeSamples cellCount = fromUnboxed1 $+ VU.replicate cellCount emptySample+++-- | Gather samples from ensemble. Return True if sampling is+-- finished, False otherwise.+updateSamples :: Ensemble+ -> MacroSamplingMonad Bool+updateSamples ens =+ let+ addCellParameters :: BasicMacroParameters+ -> BasicMacroParameters+ -> BasicMacroParameters+ addCellParameters !(n1, v1, c1) !(n2, v2, c2) =+ (n1 + n2, v1 <+> v2, c1 + c2)+ in do+ sorting@(cellCount, _) <- lift $ lift $ asks classifier++ maxSteps <- lift $ ask++ sampling <- get++ -- n is steps left for averaging+ let (n, oldSamples) =+ case sampling of+ None -> (maxSteps, initializeSamples cellCount)+ Incomplete s o -> (s, o)+ Complete _ -> error "updateSamples called, but pool's closed."+ weight = 1 / fromIntegral maxSteps+ -- Sort particles into macroscopic cells for sampling+ sorted = classifyParticles sorting ens+ -- Sampling results from current step+ stepSamples = cellMap (\_ c -> sampleMacroscopic c weight) sorted+ -- Add samples from current step to all sum of samples collected so+ -- far+ !newSamples <- R.computeP $+ R.zipWith addCellParameters oldSamples stepSamples++ let fin = (n == 0)++ -- Update state of sampling process+ put $ case fin of + True -> Complete newSamples+ False -> Incomplete (n - 1) newSamples++ return fin+++-- | Sample macroscopic values in a cell.+sampleMacroscopic :: S.Maybe CellContents+ -> Double+ -- ^ Multiply all sampled parameters by this number,+ -- which is the statistical weight of one sample.+ -- Typically this is inverse to the amount of steps+ -- used for averaging.+ -> BasicMacroParameters+sampleMacroscopic !c !weight =+ case c of+ S.Nothing -> emptySample+ S.Just ens ->+ let+ -- Particle count+ n = fromIntegral $ VU.length ens+ -- Particle averaging factor+ s = 1 / n+ -- Mean absolute velocity+ m1 = (VU.foldl' (\v0 (_, v) -> v0 <+> v) (0, 0, 0) ens) .^ s+ -- Mean square thermal velocity+ c2 = (VU.foldl' (+) 0 $+ VU.map (\(_, v) ->+ let+ thrm = (v <-> m1)+ in+ (thrm .* thrm))+ ens) * s+ in+ (n * weight, m1 .^ weight, c2 * weight)
+ src/DSMC/Motion.hs view
@@ -0,0 +1,75 @@+{-# LANGUAGE BangPatterns #-}++{-|++Collisionless motion.++-}++module DSMC.Motion+ ( motion+ )++where++import Control.Parallel.Stochastic++import qualified Data.Strict.Maybe as S++import qualified Data.Array.Repa as R++import qualified Data.Vector.Unboxed as VU++import System.Random.MWC++import DSMC.Particles+import DSMC.Surface+import DSMC.Traceables hiding (trace)+import DSMC.Util++import Control.Monad.ST++-- | Sequential action to move particles and consider particle-body+-- collisions.+reflect :: GenST s+ -> Body+ -> Time+ -> Reflector s+ -> VU.Vector Particle+ -> ST s (VU.Vector Particle)+reflect g body dt reflector ens = do+ VU.forM ens $ \pcl -> do+ -- Particle after collisionless motion+ let movedPcl = move dt pcl+ case (hitPoint dt body movedPcl) of+ -- Enjoy your convex-only case.+ S.Just (HitPoint th (S.Just n)) ->+ let+ -- Position and velocity at hit point+ (pos', v) = move th pcl+ in do+ -- Sample velocity for reflected particle+ vR <- reflector g n v+ -- Move particle away from surface with new velocity+ return $ move (-th) (pos', vR)+ _ -> return $ movedPcl+++-- | Collisionless motion step.+motion :: ParallelSeeds+ -> Body+ -> Time+ -> Surface+ -> Ensemble+ -> (Ensemble, ParallelSeeds)+motion gs b dt surf ens =+ let+ -- | Since 'reflect' is sequential, we split ensemble into N+ -- slices and process them in parallel.+ reflector = makeReflector surf+ !(v', newSeeds) = + splitParMapST (\g e -> reflect g b dt reflector e)+ (R.toUnboxed ens) gs+ in+ (fromUnboxed1 v', newSeeds)+
+ src/DSMC/Particles.hs view
@@ -0,0 +1,98 @@+{-# LANGUAGE BangPatterns #-}++{-|++Particles, ensembles, flow parameters.++-}++module DSMC.Particles+ ( -- * Particles+ Particle+ , move+ -- ** Particle ensembles+ , Ensemble+ , emptyEnsemble+ , ensembleSize+ , filterEnsemble+ , printEnsemble+ -- * Flows+ , Flow(..)+ , modelConcentration+ )++where++import qualified Data.Array.Repa as R+import qualified Data.Vector.Unboxed as VU++import DSMC.Util+import DSMC.Util.Vector+++-- | Gas particle with position and velocity.+type Particle = (Point, Vec3)+++-- | Linearly move particle for t time and update its position.+move :: Time -> Particle -> Particle+move !dt !(pos, v) = (pos <+> (v .^ dt), v)+{-# INLINE move #-}+++-- | Flow with given concentration, temperature, mass of molecule and+-- macroscopic velocity.+data Flow = Flow { concentration :: !Double+ , temperature :: !Double+ , mass :: !Double+ , velocity :: !Vec3+ , statWeight :: !Double+ -- ^ How many real particles a single simulator+ -- represents.+ }+ deriving (Show)+++-- | Calculate what model concentration will simulate real flow+-- concentration wrt statistical weight of single particle.+modelConcentration :: Flow -> Double+modelConcentration flow = (concentration flow) / (statWeight flow)+++-- | Repa array of particles.+type Ensemble = R.Array R.U R.DIM1 Particle+++-- | Ensemble with zero particles in it.+emptyEnsemble :: Ensemble+emptyEnsemble = fromUnboxed1 $ VU.empty+++-- | Amount of particles in ensemble.+ensembleSize :: Ensemble -> Int+ensembleSize ens = n where (R.Z R.:. n) = R.extent ens++-- | Print particles, one per row, using the format:+--+-- > x y z u v w+--+-- where @x y z@ are position coordinates and @u v w@ are velocity+-- components.+printEnsemble :: Ensemble -> IO ()+printEnsemble particles = do+ VU.forM_ (R.toUnboxed particles)+ (\((x, y, z), (u, v, w)) -> putStrLn $ unwords (map show [x, y, z, u, v, w]))+++-- | Filter out those particles which do not satisfy the predicate.+filterEnsemble :: (Particle -> Bool) -> Ensemble -> Ensemble+filterEnsemble pred' ens =+ let+ (R.Z R.:. size) = R.extent ens+ getter :: Int -> Particle+ getter !i = (R.!) ens (R.ix1 i)+ {-# INLINE getter #-}+ predI :: Int -> Bool+ predI !i = pred' $ getter i+ in+ R.selectP predI getter size ens
+ src/DSMC/Surface.hs view
@@ -0,0 +1,103 @@+{-# LANGUAGE BangPatterns #-}++{-|++Gas-surface interaction models.++-}++module DSMC.Surface+ ( Reflector+ , Surface(..)+ , makeReflector+ )++where++import Control.Monad.ST++import System.Random.MWC+import System.Random.MWC.Distributions (normal)++import DSMC.Util.Constants+import DSMC.Util.Vector+++-- | A function which takes PRNG state, molecular velocity, surface+-- normal and samples post-collisional wrt to impregnable wall+-- boundary condition.+type Reflector s = GenST s -> Vec3 -> Vec3 -> ST s Vec3+++-- | Surface models.+data Surface = CLL { bodyTemperature :: !Double+ -- ^ Absolute temperature of surface.+ , alpha :: !Double+ -- ^ Kinetic energy accomodation for normal+ -- velocity component.+ , sigma :: !Double+ -- ^ Accomodation for tangential momentum.+ } |+ -- ^ Cercignani-Lampis-Lord model.+ Diffuse { bodyTemperature :: !Double+ -- ^ Absolute temperature of surface.+ , mass :: !Double+ -- ^ Mass of reflected molecules (usually equal+ -- to that in incident flow).+ } |+ -- ^ Diffuse reflection.+ Mirror+ -- ^ Surface with specular reflection.+++-- | Produce reflector depending on surface type.+makeReflector :: Surface -> Reflector s+makeReflector (CLL t alphanor sigmatan) =+ let+ f = sqrt (2 * t * unigas)+ alphatan = sigmatan * (2 - sigmatan)+ cll :: Reflector s+ cll !g !n !vel =+ let+ !e1 = normalize $ n >< vel+ !e2 = normalize $ n >< e1+ !ui = -(vel .* n)+ !vi = vel .* e2+ !urm = ui * sqrt (1 - alphanor) / f+ !vrm = vi * (1 - sigmatan)+ in do+ !angle <- uniformR (0.0, pi * 2) g+ !angle' <- uniformR (0.0, pi * 2) g+ !v2 <- uniform g+ !v2' <- uniform g+ let !r = sqrt (- (alphatan * log v2))+ !vr = r * (cos angle) * f + vrm+ !wr = r * (sin angle) * f+ !r' = sqrt (- (alphanor * log v2'))+ !ur = sqrt (urm * urm + 2 * r' * urm * (cos angle') + r' * r') * f++ return $! ((e1 .^ wr) <+> (e2 .^ vr) <+> (n .^ ur))+ {-# INLINE cll #-}+ in+ cll++makeReflector Mirror = \_ !v !n -> return $! v <-> (n .^ (v .* n) .^ 2)++makeReflector (Diffuse t m) =+ let+ s = sqrt $ boltzmann * t / m+ beta = 1 / (s * (sqrt 2))+ diffuse :: Reflector s+ diffuse g n vel =+ let+ e1 = normalize $ n >< vel+ e2 = normalize $ n >< e1+ in do+ w <- normal 0 s g+ v <- normal 0 s g+ u' <- uniform g+ let u = (sqrt (- (log u'))) / beta+ return $ (e1 .^ v) <+> (e2 .^ w) <+> (n .^ u)+ {-# INLINE diffuse #-}+ in+ diffuse
+ src/DSMC/Traceables.hs view
@@ -0,0 +1,486 @@+{-# LANGUAGE BangPatterns #-}++{-|++Ray-casting routines for constructive solid geometry.++This module provides constructors for complex bodies as well as+routines to compute intersections of such bodies with ray. In DSMC it+is used to calculate points at which particles hit the body surface.++Gas-surface interactions are not handled by this module, see+'DSMC.Surface' instead.++-}++module DSMC.Traceables+ ( -- * Bodies+ Body+ -- ** Primitives+ , plane+ , sphere+ , cylinder+ , cylinderFrustum+ , cone+ , coneFrustum+ -- ** Compositions+ , intersect+ , unite+ , complement+ -- * Ray casting+ , HitPoint(..)+ , hitPoint+ , HitSegment+ , Trace+ , trace+ -- * Body membership+ , inside+ )++where++import Prelude hiding (Just, Nothing, Maybe, fst, snd)++import Data.Functor+import Data.Strict.Maybe+import Data.Strict.Tuple++import DSMC.Particles+import DSMC.Util+import DSMC.Util.Vector+++-- | Time when particle hits the surface with normal at the hit point.+-- If hit is in infinity, then normal is Nothing.+--+-- Note that this datatype is strict only on first argument: we do not+-- compare normals when classifying traces.+data HitPoint = HitPoint !Double (Maybe Vec3)+ deriving (Eq, Ord, Show)+++-- | A segment on time line when particle is inside the body.+--+-- Using strict tuple performs better: 100 traces for 350K+-- particles perform roughly 7s against 8s with common datatypes.+type HitSegment = (Pair HitPoint HitPoint)+++-- | Trace of a linearly-moving particle on a body is a list of time+-- segments/intervals during which the particle is inside the body.+--+-- > # - particle+-- > \+-- > \+-- > o------------+-- > ---/ * \---+-- > -/ * \-+-- > / * \+-- > ( * - trace )+-- > \ * /+-- > -\ * /-+-- > primitive - ---\ * /---+-- > --------o----+-- > \+-- > \+-- > _\/+-- > \+--+--+-- For example, since a ray intersects a plane only once, a half-space+-- primitive defined by this plane results in a half-interval trace of+-- a particle:+--+-- > /+-- > /+-- > /+-- > #------------------------o*****************>+-- > | / |+-- > particle / goes to infinity+-- > /+-- > /+-- > /+-- > / - surface of half-space+--+-- Ends of segments or intervals are calculated by intersecting the+-- trajectory ray of a particle and the surface of the primitive. This+-- may be done by substituting the equation of trajectory @X(t) = X_o+-- + V*t@ into the equation which defines the surface and solving it+-- for @t@. If the body is a composition, traces from primitives are+-- then classified according to operators used to define the body+-- (union, intersection or complement).+--+-- Although only convex primitives are used in current implementation,+-- compositions may result in concave bodies, which is why trace is+-- defined as a list of segments.+--+--+-- In this example, body is an intersection of a sphere and sphere+-- complement:+--+-- > /|\+-- > |+-- > |+-- > |+-- > ----------- |+-- > ----/ \--o-+-- > -/ * \-+-- > -/ hs2 - * \+-- > -/ * ---/+-- > / o/+-- > / -/|+-- > / / |+-- > | / |+-- > / | |+-- > | / |+-- > | | |+-- > | \ |+-- > \ | |+-- > | \ |+-- > \ \ |+-- > \ -\|+-- > \ o\+-- > -\ * ---\+-- > -\ hs1 - * /+-- > -\ * /-+-- > ----\ /--o-+-- > ----------- |+-- > |+-- > |+-- > # - particle+--+-- If only intersections of concave primitives were allowed, then+-- trace type might be simplified to be just single 'HitSegment'.+type Trace = [HitSegment]+++-- | IEEE positive infinity.+infinityP :: Double+infinityP = (/) 1 0+++-- | Negative infinity.+infinityN :: Double+infinityN = -infinityP+++hitN :: HitPoint+hitN = (HitPoint infinityN Nothing)+++hitP :: HitPoint+hitP = (HitPoint infinityP Nothing)+++-- | CSG body is a recursive composition of primitive objects or other+-- bodies.+data Body = Plane !Vec3 !Double+ -- ^ Half-space with normalized outward normal and distance+ -- of boundary plane from origin.+ | Sphere !Vec3 !Double+ -- ^ Sphere defined by center and radius.+ | Cylinder !Vec3 !Point !Double+ -- ^ Infinite circular cylinder with normalized axis vector,+ -- point on axis and radius.+ | Cone !Vec3 !Point !Double !Matrix !Double !Double+ -- ^ Cone defined by inward axis direction, vertex and+ -- cosine to angle h between axis and outer edge.+ --+ -- Additionally transformation matrix $n * n - cos^2 h$,+ -- tangent of angle and odelta are stored for intersection+ -- calculations.+ | Union !Body !Body+ | Intersection !Body !Body+ | Complement !Body+ deriving Show+++-- | A half-space defined by arbitary point on the boundary plane and+-- outward normal (not necessarily a unit vector).+plane :: Point -> Vec3 -> Body+plane p n = Plane nn (p .* nn)+ where+ nn = normalize n+++-- | A sphere defined by center point and radius.+sphere :: Vec3 -> Double -> Body+sphere o r = Sphere o r+++-- | An infinite circular cylinder defined by two arbitary+-- points on axis and radius.+cylinder :: Point -> Point -> Double -> Body+cylinder p1 p2 r = Cylinder (normalize $ p2 <-> p1) p1 r+++-- | A finite right circular cylinder defined by two points on its top+-- and bottom and radius.+cylinderFrustum :: Point -> Point -> Double -> Body+cylinderFrustum pb pt r =+ intersect (plane pt axis)+ (intersect (plane pb $ invert axis)+ (cylinder pb pt r))+ where+ axis = pt <-> pb+++-- | An infinite right circular cone defined by outward axis vector,+-- apex point and angle between generatrix and axis (in degrees, less+-- than 90).+cone :: Vec3 -> Point -> Double -> Body+cone a o h =+ let+ h' = cos $ (h * pi / 180)+ n = normalize $ invert a+ gamma = diag (-h' * h')+ m = addM (n `vxv` n) gamma+ ta = tan $ h+ odelta = n .* o+ in+ Cone n o h' m ta odelta+++-- | A conical frustum given by two points on its axis with radii at+-- that points. One of radii may be zero (in which case one of frustum+-- ends will be the apex).+coneFrustum :: (Point, Double) -> (Point, Double) -> Body+coneFrustum (p1, r1) (p2, r2) =+ let+ -- Direction from pb to pt is towards apex. Corresponding+ -- radii are rb > rt.+ (pb, rb, pt, rt) = case (r1 < r2) of+ True -> (p2, r2, p1, r1)+ False -> (p1, r1, p2, r2)+ -- Cone axis and frustum height+ gap = pt <-> pb+ height = norm gap+ axis = normalize gap+ -- Calculate distance from pt to apex.+ dist = if rt == 0 + then 0 + else height / (rb / rt - 1)+ apex = pt <+> (axis .^ dist)+ -- Angle between generatrix and axis+ degs = atan (rb / (dist + norm (pt <-> pb))) * (180 / pi)+ in+ intersect (plane pt axis)+ (intersect (plane pb $ invert axis)+ (cone axis apex degs))+++-- | Intersection of two bodies.+intersect :: Body -> Body -> Body+intersect !b1 !b2 = Intersection b1 b2+++-- | Union of two bodies.+unite :: Body -> Body -> Body+unite !b1 !b2 = Union b1 b2+++-- | Complement to a body (normals flipped).+complement :: Body -> Body+complement !b = Complement b+++-- | Calculate a trace of a particle on a body.+trace :: Body -> Particle -> Trace+{-# INLINE trace #-}++trace !b@(Plane n d) !p@(pos, v) =+ let+ !f = -(n .* v)+ in+ if f == 0+ then+ -- If ray is parallel to plane and is inside, then trace is+ -- the whole timeline.+ if inside b p+ then [(HitPoint infinityN Nothing) :!: (HitPoint infinityP Nothing)]+ else []+ else+ let+ !t = (pos .* n - d) / f+ in+ if f > 0+ then [(HitPoint t (Just n)) :!: (HitPoint infinityP Nothing)]+ else [(HitPoint infinityN Nothing) :!: (HitPoint t (Just n))]++trace !(Sphere c r) !(pos, v) =+ let+ !d = pos <-> c+ !roots = solveq (v .* v) (v .* d * 2) (d .* d - r * r)+ normal !u = normalize (u <-> c)+ in+ case roots of+ Nothing -> []+ Just (t1 :!: t2) ->+ [HitPoint t1 (Just $ normal $ moveBy pos v t1) :!:+ HitPoint t2 (Just $ normal $ moveBy pos v t2)]++trace !(Cylinder n c r) !(pos, v) =+ let+ d = (pos <-> c) >< n+ e = v >< n+ roots = solveq (e .* e) (d .* e * 2) (d .* d - r * r)+ normal u = normalize $ h <-> (n .^ (h .* n))+ where h = u <-> c+ in+ case roots of+ Nothing -> []+ Just (t1 :!: t2) ->+ [HitPoint t1 (Just $ normal $ moveBy pos v t1) :!:+ HitPoint t2 (Just $ normal $ moveBy pos v t2)]++trace !(Cone n c _ m ta odelta) !(pos, v) =+ let+ delta = pos <-> c+ c2 = dotM v v m+ c1 = dotM v delta m+ c0 = dotM delta delta m+ roots = solveq c2 (2 * c1) c0+ normal !u = normalize $ nx .^ (1 / ta) <-> ny .^ ta+ where h = u <-> c+ -- Component of h parallel to cone axis+ ny' = n .^ (n .* h)+ ny = normalize ny'+ -- Perpendicular component+ nx = normalize $ h <-> ny'+ in+ case roots of+ Nothing -> []+ Just (t1 :!: t2) ->+ let+ pos1 = moveBy pos v t1+ pos2 = moveBy pos v t2+ in+ case ((pos1 .* n - odelta) > 0, (pos2 .* n - odelta) > 0) of+ (True, True) -> [HitPoint t1 (Just $ normal pos1) :!:+ HitPoint t2 (Just $ normal pos2)]+ (True, False) -> [HitPoint infinityN Nothing :!:+ HitPoint t1 (Just $ normal pos1)]+ (False, True) -> [HitPoint t2 (Just $ normal pos2) :!:+ HitPoint infinityP Nothing]+ (False, False) -> []++trace !(Intersection b1 b2) !p =+ intersectTraces tr1 tr2+ where+ tr1 = trace b1 p+ tr2 = trace b2 p++trace !(Union b1 b2) !p =+ uniteTraces tr1 tr2+ where+ tr1 = trace b1 p+ tr2 = trace b2 p++trace !(Complement b) !p =+ complementTrace $ trace b p+++uniteTraces :: Trace -> Trace -> Trace+uniteTraces u [] = u+uniteTraces u (v:t2) =+ uniteTraces (unite1 u v) t2+ where+ merge :: HitSegment -> HitSegment -> HitSegment+ merge !(a1 :!: b1) !(a2 :!: b2) = (min a1 a2) :!: (max b1 b2)+ {-# INLINE merge #-}+ unite1 :: Trace -> HitSegment -> Trace+ unite1 [] hs = [hs]+ unite1 t@(hs1@(a1 :!: b1):tr') hs2@(a2 :!: b2)+ | b1 < a2 = hs1:(unite1 tr' hs2)+ | a1 > b2 = hs2:t+ | otherwise = unite1 tr' (merge hs1 hs2)+ {-# INLINE unite1 #-}+{-# INLINE uniteTraces #-}+++intersectTraces :: Trace -> Trace -> Trace+intersectTraces tr1 tr2 =+ let+ -- Overlap two overlapping segments+ overlap :: HitSegment -> HitSegment -> HitSegment+ overlap !(a1 :!: b1) !(a2 :!: b2) = (max a1 a2) :!: (min b1 b2)+ {-# INLINE overlap #-}+ in+ case tr2 of+ [] -> []+ (hs2@(a2 :!: b2):tr2') ->+ case tr1 of+ [] -> []+ (hs1@(a1 :!: b1):tr1') ->+ case (b1 < a2) of+ True -> (intersectTraces tr1' tr2)+ False ->+ case (b2 < a1) of+ True -> intersectTraces tr1 tr2'+ False -> (overlap hs1 hs2):(intersectTraces tr1' tr2)+{-# INLINE intersectTraces #-}+++-- | Complement to trace (normals flipped) in @R^3@.+complementTrace :: Trace -> Trace+complementTrace ((sp@(HitPoint ts _) :!: ep):xs) =+ start ++ (complementTrace' ep xs)+ where+ flipNormals :: HitSegment -> HitSegment+ flipNormals !((HitPoint t1 n1) :!: (HitPoint t2 n2)) =+ (HitPoint t1 (invert <$> n1)) :!: (HitPoint t2 (invert <$> n2))+ {-# INLINE flipNormals #-}+ -- Start from infinity if first hitpoint is finite+ start = if (isInfinite ts)+ then []+ else [flipNormals $ hitN :!: sp]+ complementTrace' :: HitPoint -> Trace -> Trace+ complementTrace' c ((a :!: b):tr) =+ -- Bridge between last point of previous segment and first+ -- point of the next one.+ (flipNormals (c :!: a)):(complementTrace' b tr)+ complementTrace' a@(HitPoint t _) [] =+ -- End in infinity if last hitpoint is finite+ if (isInfinite t)+ then []+ else [flipNormals (a :!: hitP)]+complementTrace [] = [hitN :!: hitP]+{-# INLINE complementTrace #-}+++-- | If the particle has hit the body during last time step, calculate+-- the first corresponding 'HitPoint'. Note that the time at which the hit+-- occured will be negative. This is the primary function to calculate+-- ray-body intersections.+hitPoint :: Time -> Body -> Particle -> Maybe HitPoint+hitPoint !dt !b !p =+ let+ lastHit = [(HitPoint (-dt) Nothing) :!: (HitPoint 0 Nothing)]+ in+ case (intersectTraces lastHit $ trace b p) of+ [] -> Nothing+ (hs:_) -> Just $ fst hs+{-# INLINE hitPoint #-}+++-- | True if particle is in inside the body.+inside :: Body -> Particle -> Bool+{-# INLINE inside #-}++inside !(Plane n d) !(pos, _) = (pos .* n - d) < 0++inside !(Sphere c r) !(pos, _) = (norm $ pos <-> c) < r++inside !(Cylinder n c r) !(pos, _) =+ (norm $ h <-> (n .^ (h .* n))) < r+ where+ h = pos <-> c++inside !(Cone n c a _ _ _) !(pos, _) =+ (n .* (normalize $ pos <-> c)) > a++inside !(Intersection b1 b2) !p = inside b1 p && inside b2 p++inside !(Union b1 b2) !p = inside b1 p || inside b2 p++inside !(Complement b) !p = not $ inside b p
+ src/DSMC/Traceables/Parser.hs view
@@ -0,0 +1,288 @@+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE ScopedTypeVariables #-}++-- | Parser for body definitions.+--+-- Body definition contains a number of solid definitions and ends+-- with the top level object definition. RHS of solid equations may+-- reference other solids to compose into complex bodies.+--+-- Multiple-body compositions are right-associative.+--+-- > # comment+-- >+-- > # define few primitives+-- > solid b1 = sphere (0, 0, 0; 5);+-- > solid p1 = plane (0, 0, 0; 1, 0, 0);+-- >+-- > # define a composition+-- > solid body = b1 and p1;+-- >+-- > # assign it to be the top level object+-- > tlo body;+--+-- Statements must end with a semicolon (newlines are optional).+-- Excessive spaces are ignored.+--+-- Top-level object line must reference a previously defined solid.+--+-- Syntax for primitives follows the signatures of 'Traceables'+-- constructors for 'T.plane' and 'T.sphere', but differs for cylinder+-- and cone, as this module provides access only to frustums+-- ('T.cylinderFrustum' and 'T.coneFrustum').+--+-- [Half-space] @plane (px, py, pz; nx, ny, nz)@, where @(px, py, pz)@+-- is a point on a plane which defines the half-space and @(nx, ny,+-- nz)@ is a normal to the plane (outward to the half-space), not+-- necessarily a unit vector.+--+-- [Sphere] @sphere (cx, cy, cz; r)@, where @(cx, cy, cz)@ is a+-- central point of a sphere and @r@ is radius.+--+-- [Right circular cylinder] @cylinder (p1x, p1y, p1z; p2x, p2y, p2z;+-- r)@ where @(p1x, p1y, p1z)@ and @(p2x, p2y, p2z)@ are bottom and+-- top points on axis and @r@ is radius.+--+-- [Right circular conical frustum] @cone (p1x, p1y, p1z; r1; p2x,+-- p2y, p2z; r2)@ where @(p1x, p1y, p1z)@ and @(p2x, p2y, p2z)@ are+-- bottom and top points on cone axis and @r1@, @r2@ are the+-- corresponding radii.++module DSMC.Traceables.Parser+ ( parseBody+ , parseBodyFile+ )++where++import Prelude as P++import Control.Applicative+import qualified Control.Exception as E+import Control.Monad++import Control.Monad.Trans.Class+import Control.Monad.Trans.State.Strict++import Data.Attoparsec.Char8+import Data.ByteString.Char8 as B++import qualified Data.Map as M++import qualified DSMC.Traceables as T+import DSMC.Util.Vector+++-- | Transformer which adds lookup table to underlying monad.+type Table a k v = StateT (M.Map k v) a+++-- | Add entry to the lookup table.+addEntry :: (Ord k, Monad a) => k -> v -> Table a k v ()+addEntry key value = liftM (M.insert key value) get >>= put+++-- | Lookup entry in the table.+getEntry :: (Ord k, Monad a) => k -> Table a k v (Maybe v)+getEntry key = liftM (M.lookup key) get+++-- | Parser with lookup table.+type CSGParser = Table Parser String T.Body+++lp :: Parser Char+lp = char '('+++rp :: Parser Char+rp = char ')'+++eq :: Parser Char+eq = char '='+++cancer :: Parser Char+cancer = char ';'+++comma :: Parser Char+comma = char ','+++-- | Read comma-separated three doubles into point.+--+-- > <triple> ::= <double> ',' <double> ',' <double>+triple :: Parser Point+triple = (,,) <$> double+ <*>+ (skipSpace *> comma *> skipSpace *>+ double+ <* skipSpace <* comma <* skipSpace)+ <*>+ double+++keywords :: [String]+keywords = [ "solid"+ , "tlo"+ , "plane"+ , "sphere"+ , "cylinder"+ , "cone"+ ]+++-- | Read variable name or fail if it's a keyword.+varName :: CSGParser String+varName = do+ k <- lift $ many1 (letter_ascii <|> digit)+ case (P.elem k keywords) of+ False -> return k+ True -> fail $ "Unexpected keyword: " ++ k+++-- | Lookup body in table by its name or fail if it is undefined.+readName :: CSGParser T.Body+readName = do+ k <- varName+ v <- getEntry k+ case v of+ Just b -> return b+ _ -> fail $ "Undefined solid: " ++ k+++-- > <plane> ::=+-- > 'plane (' <triple> ';' <triple> ')'+plane :: Parser T.Body+plane = T.plane <$>+ (string "plane" *> skipSpace *> lp *> skipSpace *> triple) <*>+ (skipSpace *> cancer *> skipSpace *> triple <* skipSpace <* rp)+++-- > <sphere> ::=+-- > 'sphere (' <triple> ';' <double> ')'+sphere :: Parser T.Body+sphere = T.sphere <$>+ (string "sphere" *> skipSpace *> lp *> skipSpace *> triple) <*>+ (skipSpace *> cancer *> skipSpace *> double <* skipSpace <* rp)+++-- > <cylinder> ::=+-- > 'cylinder (' <triple> ';' <triple> ';' <double> ')'+cylinder :: Parser T.Body+cylinder = T.cylinderFrustum <$>+ (string "cylinder" *> skipSpace *> lp *> skipSpace *> triple) <*>+ (skipSpace *> cancer *> skipSpace *> triple) <*>+ (skipSpace *> cancer *> skipSpace *> double <* skipSpace <* rp)+++-- > <cone> ::=+-- > 'cone (' <triple> ';' <double> ';' <triple> ';' <double> ')'+cone :: Parser T.Body+cone = T.coneFrustum <$> + ((,) <$>+ (string "cone" *> skipSpace *> lp *> skipSpace *> triple) <*>+ (skipSpace *> cancer *> skipSpace *> double)) <*>+ ((,) <$>+ (skipSpace *> cancer *> skipSpace *> triple) <*>+ (skipSpace *> cancer *> skipSpace *> double <* skipSpace <* rp))+++primitive :: Parser T.Body+primitive = plane <|> sphere <|> cylinder <|> cone+++-- > <complement> ::= 'not' <body>+complement :: CSGParser T.Body+complement = T.complement <$> (lift (string "not" *> skipSpace) *> body)+++-- > <union> ::= <uncomposed-body> 'or' <body>+union :: CSGParser T.Body+union = binary "or" T.unite+++-- > <intersection> ::= <uncomposed-body> 'and' <body>+intersection :: CSGParser T.Body+intersection = binary "and" T.intersect+++-- | Parse binary operation on two bodies with given composition+-- operators.+binary :: ByteString -> (T.Body -> T.Body -> T.Body) -> CSGParser T.Body+binary op compose = do+ b1 <- uncomposedBody+ lift (skipSpace *> string op *> skipSpace)+ b2 <- body+ return $ compose b1 b2+++-- | Read stamement which adds new solid entry to lookup table.+--+-- > <statement> ::=+-- > 'solid' <varname> '=' <body> ';'+statement :: CSGParser ()+statement = do+ lift $ string "solid" *> skipSpace+ k <- varName+ lift $ skipSpace <* eq <* skipSpace+ v <- body <* lift (cancer *> skipSpace)+ addEntry k v+++-- | Expression is either a primitive, a reference to previously+-- defined solid or an operation on expressions.+--+-- > <body> ::= <union> | <intersection> | <complement> | <primitive> | <reference>+body :: CSGParser T.Body+body = union <|> intersection <|> complement <|> uncomposedBody+++-- Used to terminate left branch of binary compositions.+--+-- > <uncomposed-body> ::= <primitive> | <reference>+uncomposedBody :: CSGParser T.Body+uncomposedBody = lift primitive <|> readName+++-- | Top-level object declaration.+--+-- > <tlo> ::= 'tlo' <body> ';'+topLevel :: CSGParser T.Body+topLevel = lift (string "tlo" *> skipSpace) *>+ readName+ <* lift (cancer <* skipSpace)+++-- | Read one-line comment starting with hash sign.+comment :: Parser ()+comment = char '#' >> (manyTill anyChar endOfLine) >> return ()+++-- | Read sequence of statements which define solids, and finally read+-- top level object definition.+--+-- > <geoFile> ::= <statement> <geoFile> | <comment> <geoFile> | <tlo>+geoFile :: CSGParser T.Body+geoFile = (many1 $ lift comment <|> statement) *> topLevel+++-- | Try to read body definition from bytestring. Return body or error+-- message if parsing fails.+parseBody :: ByteString -> Either String T.Body+parseBody input =+ case (parseOnly (runStateT geoFile M.empty) input) of+ Right (b, _) -> Right b+ Left msg -> Left msg+++-- | Read body definition from file. If parsing fails or IOError when+-- reading file occurs, return error message.+parseBodyFile :: FilePath -> IO (Either String T.Body)+parseBodyFile file = do+ res <- E.try $ B.readFile file+ return $ case res of+ Right d -> parseBody d+ Left e -> Left $ show (e :: E.IOException)
+ src/DSMC/Util.hs view
@@ -0,0 +1,107 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE Rank2Types #-}++-- | Utility functions and definitions.++module DSMC.Util+ ( solveq+ , SquareRoots+ , fromUnboxed1+ , iforM_+ , randomSeed+ , Time+ , DSMCRootMonad+ )++where++import Prelude hiding (Just, Nothing, Maybe, fst)++import Data.Bits+import Data.Word++import Data.Functor+import Data.List++import qualified Data.ByteString as BS++import Data.Strict.Maybe+import Data.Strict.Tuple++import qualified Data.Array.Repa as R+import qualified Data.Vector.Unboxed as VU++import System.Entropy+import System.Random.MWC++import Data.Splittable++-- | Results of solving a quadratic equation.+type SquareRoots = Maybe (Pair Double Double)++-- | Solve quadratic equation @ax^2 + bx + c = 0@.+--+-- If less than two roots exist, Nothing is returned.+solveq :: Double+ -- ^ a+ -> Double+ -- ^ b+ -> Double+ -- ^ c+ -> SquareRoots+solveq !a !b !c+ | (d > 0) = Just $ min r1 r2 :!: max r1 r2+ | otherwise = Nothing+ where+ d = b * b - 4 * a * c+ q = sqrt d+ t = 2 * a+ r = - b / t+ s = q / t+ r1 = r - s+ r2 = r + s+{-# INLINE solveq #-}+++-- | Convert between Repa 'R.DIM1'-arrays and unboxed 'VU.Vector's.+fromUnboxed1 :: (VU.Unbox e) => VU.Vector e -> R.Array R.U R.DIM1 e+fromUnboxed1 v = R.fromUnboxed (R.ix1 $ VU.length v) v+{-# INLINABLE fromUnboxed1 #-}+++-- | Map monadic action over pairs of vector indices and items and+-- throw away the results.+iforM_ :: (Monad m, VU.Unbox a) =>+ VU.Vector a+ -> ((Int, a) -> m b)+ -> m ()+iforM_ v = VU.forM_ (VU.imap (,) v)+{-# INLINE iforM_ #-}+++-- | Fetch vector of random Word32 values from system entropy source.+randomWord32Vector :: Int -> IO (VU.Vector Word32)+randomWord32Vector n =+ let+ -- Left fold accumulator to shift Word8 onto Word32+ accum a o = (a `shiftL` 8) .|. fromIntegral o+ in do+ w8stream <- getEntropy (n * 4)+ -- Split the stream into 4-length lists of Word8+ let w8s = map BS.unpack $ splitIn n w8stream+ return $ VU.fromList $ map (foldl' accum 0) w8s+++-- | Fetch new RNG seed from system entropy source.+randomSeed :: IO Seed+randomSeed = toSeed <$> randomWord32Vector 256+++-- | Time in seconds.+type Time = Double+++-- | Several modules define a chain of monads to maintain context of+-- the running simulation. In its root is the IO monad which we use to+-- send logger messages from monads atop the root.+type DSMCRootMonad = IO
+ src/DSMC/Util/Constants.hs view
@@ -0,0 +1,31 @@+-- | Physical constants.++module DSMC.Util.Constants+ ( amu+ , avogadro+ , boltzmann+ , unigas+ )++where+++-- | Atomic mass unit 1.660538921(73)e-27, inverse to Avogadro's+-- constant.+amu :: Double+amu = 1.6605389217373737e-27+++-- | Avogadro constant 6.02214129(27)e23+avogadro :: Double+avogadro = 6.0221412927272727e23+++-- | Boltzmann constant 1.3806488(13)e-23+boltzmann :: Double+boltzmann = 1.3806488131313131e-23+++-- | Universal gas constant.+unigas :: Double+unigas = boltzmann * avogadro
+ src/DSMC/Util/Vector.hs view
@@ -0,0 +1,169 @@+{-# LANGUAGE BangPatterns #-}++{-| ++Simple 3-vectors and matrices built atop tuples.++-}++module DSMC.Util.Vector+ ( Vec3+ , Matrix+ , Point+ , origin+ -- * Vector operations+ , (<+>)+ , (<->)+ , (><)+ , (.^)+ , (.*)+ , norm+ , normalize+ , invert+ , distance+ , moveBy+ -- * Matrix operations+ , mxv+ , vxv+ , dotM+ , diag+ , addM+ -- * Cartesian system+ , buildCartesian+ )++where++import Prelude hiding (reverse)+++-- | Vector in @R^3@.+type Vec3 = (Double, Double, Double)+++-- | Matrix given by its rows.+type Matrix = (Vec3, Vec3, Vec3)+++-- | Point in @R^3@.+type Point = Vec3+++-- | Origin point @(0, 0, 0)@.+origin :: Point+origin = (0, 0, 0)+++-- | Add two vectors.+(<+>) :: Vec3 -> Vec3 -> Vec3+(<+>) !(x1, y1, z1) !(x2, y2, z2) = (x1 + x2, y1 + y2, z1 + z2)+{-# INLINE (<+>) #-}+++-- | Subtract two vectors.+(<->) :: Vec3 -> Vec3 -> Vec3+(<->) !(x1, y1, z1) !(x2, y2, z2) = (x1 - x2, y1 - y2, z1 - z2)+{-# INLINE (<->) #-}+++-- | Vec3 cross product.+(><) :: Vec3 -> Vec3 -> Vec3+(><) !(x1, y1, z1) !(x2, y2, z2) =+ (y1 * z2 - y2 * z1, x2 * z1 - x1 * z2, x1 * y2 - x2 * y1)+{-# INLINE (><) #-}+++-- | Scale vector.+(.^) :: Vec3 -> Double -> Vec3+(.^) !(x, y, z) !s = (x * s, y * s, z * s)+{-# INLINE (.^) #-}+++-- | Vec3 dot product.+(.*) :: Vec3 -> Vec3 -> Double+(.*) !(x1, y1, z1) !(x2, y2, z2) = x1 * x2 + y1 * y2 + z1 * z2+{-# INLINE (.*) #-}+++-- | Generic vector dot product.+--+-- Multiply transpose of first vector by given matrix, then multiply+-- the result by second vector.+dotM :: Vec3 -> Vec3 -> Matrix -> Double+dotM !v1 !v2 !m = v1 .* (m `mxv` v2)+{-# INLINE dotM #-}+++-- | Multiply matrix (given by row vectors) and vector+mxv :: Matrix -> Vec3 -> Vec3+mxv !(r1, r2, r3) !v = (r1 .* v, r2 .* v, r3 .* v)+{-# INLINE mxv #-}+++-- | Produce matrix with diagonal elements equal to given value.+diag :: Double -> Matrix+diag !d = ((d, 0, 0), (0, d, 0), (0, 0, d))+{-# INLINE diag #-}+++-- | Transpose vector and multiply it by another vector, producing a+-- matrix.+vxv :: Vec3 -> Vec3 -> Matrix+vxv !(v11, v12, v13) !v2 = (v2 .^ v11, v2 .^ v12, v2 .^ v13)+{-# INLINE vxv #-}+++-- | Euclidean distance between two points.+distance :: Point -> Point -> Double+distance !v1 !v2 = norm (v1 <-> v2)+{-# INLINE distance #-}+++-- | Euclidean norm of vector.+norm :: Vec3 -> Double+norm !(x, y, z) = sqrt (x * x + y * y + z * z)+{-# INLINE norm #-}+++-- | Produce unit vector with same direction as the original one.+normalize :: Vec3 -> Vec3+normalize !v = v .^ (1 / norm v)+{-# INLINE normalize #-}+++-- | Scale vector by -1.+invert :: Vec3 -> Vec3+invert !v = v .^ (-1)+{-# INLINE invert #-}+++-- | Move point by velocity vector for given time and return new+-- position.+moveBy :: Point+ -- ^ Current position.+ -> Vec3+ -- ^ Velocity.+ -> Double + -- ^ Time step.+ -> Point+moveBy !p !v !t = p <+> (v .^ t)+{-# INLINE moveBy #-}+++-- | Add two matrices.+--+-- We could add Applicative instance for Matrix and lift (+) to it.+addM :: Matrix -> Matrix -> Matrix+addM !(r11, r12, r13) !(r21, r22, r23) =+ (r11 <+> r21, r12 <+> r22, r13 <+> r23)+{-# INLINE addM #-}+++-- | Build cartesian axes from yaw and pitch with 0 roll. Angles are+-- in radians.+buildCartesian :: Double -> Double -> (Vec3, Vec3, Vec3)+buildCartesian yaw pitch = (u, v, w)+ where u = (cos yaw * cos pitch, sin yaw * cos pitch, sin pitch)+ v = (- (sin yaw), cos yaw, 0)+ w = u >< v+{-# INLINE buildCartesian #-}
+ src/Data/Splittable.hs view
@@ -0,0 +1,78 @@+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE UndecidableInstances #-}++module Data.Splittable + ( Split(..)+ , Combine(..)+ )+ +where++import Prelude hiding (splitAt)+import qualified Prelude as P++import qualified Data.ByteString as BS++import qualified Data.Vector.Generic as VG++ +-- | Class of tasks (containers which hold input data) which may be+-- splitted into subtasks.+--+-- 'splitAt' and 'splitIn' must preserve linear ordering on elements+-- of task, if there's any.+class Split source where+ -- | Split the source in two subsources given the size of the first source.+ splitAt :: Int+ -> source+ -> (source, source)++ -- | Calculate the overall size of the source+ size :: source -> Int++ splitIn :: Int+ -- ^ Split the source into that many subsources of equal+ -- size. This number is capped to the size of source if it+ -- exceeds it.+ -> source+ -> [source]+ splitIn n l | n < 1 = error "Can't split in less than one chunk!"+ | n > (size l) = splitIn (size l) l+ | otherwise =+ let+ partSize = (size l) `div` n+ splitIn1 acc 1 rest = acc ++ [rest]+ splitIn1 acc m rest = splitIn1 (acc ++ [v1]) (m - 1) v2+ where+ (v1, v2) = splitAt partSize rest+ in+ splitIn1 [] n l+ {-# INLINE splitIn #-}+++-- | A counterpart for 'Split'. Class of containers which may be+-- combined into single container (which holds output data).+class Combine result where+ -- | Combine list of results, preserving linear ordering.+ combine :: [result] -> result+++instance Split [s] where+ splitAt = P.splitAt+ size = P.length++instance Combine [s] where+ combine = concat+++instance (VG.Vector v e) => Split (v e) where+ splitAt = VG.splitAt+ size = VG.length++instance (VG.Vector v e) => Combine (v e) where+ combine = VG.concat+++instance Split BS.ByteString where+ splitAt = BS.splitAt+ size = BS.length