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/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/dsmc.cabal b/dsmc.cabal new file mode 100644 --- /dev/null +++ b/dsmc.cabal @@ -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.* diff --git a/src/Control/Parallel/Stochastic.hs b/src/Control/Parallel/Stochastic.hs new file mode 100644 --- /dev/null +++ b/src/Control/Parallel/Stochastic.hs @@ -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] diff --git a/src/DSMC.hs b/src/DSMC.hs new file mode 100644 --- /dev/null +++ b/src/DSMC.hs @@ -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 diff --git a/src/DSMC/Cells.hs b/src/DSMC/Cells.hs new file mode 100644 --- /dev/null +++ b/src/DSMC/Cells.hs @@ -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) diff --git a/src/DSMC/Domain.hs b/src/DSMC/Domain.hs new file mode 100644 --- /dev/null +++ b/src/DSMC/Domain.hs @@ -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) diff --git a/src/DSMC/Macroscopic.hs b/src/DSMC/Macroscopic.hs new file mode 100644 --- /dev/null +++ b/src/DSMC/Macroscopic.hs @@ -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) diff --git a/src/DSMC/Motion.hs b/src/DSMC/Motion.hs new file mode 100644 --- /dev/null +++ b/src/DSMC/Motion.hs @@ -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) + diff --git a/src/DSMC/Particles.hs b/src/DSMC/Particles.hs new file mode 100644 --- /dev/null +++ b/src/DSMC/Particles.hs @@ -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 diff --git a/src/DSMC/Surface.hs b/src/DSMC/Surface.hs new file mode 100644 --- /dev/null +++ b/src/DSMC/Surface.hs @@ -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 diff --git a/src/DSMC/Traceables.hs b/src/DSMC/Traceables.hs new file mode 100644 --- /dev/null +++ b/src/DSMC/Traceables.hs @@ -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 diff --git a/src/DSMC/Traceables/Parser.hs b/src/DSMC/Traceables/Parser.hs new file mode 100644 --- /dev/null +++ b/src/DSMC/Traceables/Parser.hs @@ -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 :: 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 :: Parser T.Body +plane = T.plane <$> + (string "plane" *> skipSpace *> lp *> skipSpace *> triple) <*> + (skipSpace *> cancer *> skipSpace *> triple <* skipSpace <* rp) + + +-- > ::= +-- > 'sphere (' ';' ')' +sphere :: Parser T.Body +sphere = T.sphere <$> + (string "sphere" *> skipSpace *> lp *> skipSpace *> triple) <*> + (skipSpace *> cancer *> skipSpace *> double <* skipSpace <* rp) + + +-- > ::= +-- > 'cylinder (' ';' ';' ')' +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 :: 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 + + +-- > ::= 'not' +complement :: CSGParser T.Body +complement = T.complement <$> (lift (string "not" *> skipSpace) *> body) + + +-- > ::= 'or' +union :: CSGParser T.Body +union = binary "or" T.unite + + +-- > ::= 'and' +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. +-- +-- > ::= +-- > 'solid' '=' ';' +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 :: CSGParser T.Body +body = union <|> intersection <|> complement <|> uncomposedBody + + +-- Used to terminate left branch of binary compositions. +-- +-- > ::= | +uncomposedBody :: CSGParser T.Body +uncomposedBody = lift primitive <|> readName + + +-- | Top-level object declaration. +-- +-- > ::= 'tlo' ';' +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 :: 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) diff --git a/src/DSMC/Util.hs b/src/DSMC/Util.hs new file mode 100644 --- /dev/null +++ b/src/DSMC/Util.hs @@ -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 diff --git a/src/DSMC/Util/Constants.hs b/src/DSMC/Util/Constants.hs new file mode 100644 --- /dev/null +++ b/src/DSMC/Util/Constants.hs @@ -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 diff --git a/src/DSMC/Util/Vector.hs b/src/DSMC/Util/Vector.hs new file mode 100644 --- /dev/null +++ b/src/DSMC/Util/Vector.hs @@ -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 #-} diff --git a/src/Data/Splittable.hs b/src/Data/Splittable.hs new file mode 100644 --- /dev/null +++ b/src/Data/Splittable.hs @@ -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