module DSMC.Domain
( Domain(..)
, getDimensions
, getCenter
, makeDomain
, initializeParticles
, openBoundaryInjection
, DomainSeeds
, clipToDomain
, 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
data Domain = Domain !Double !Double !Double !Double !Double !Double
deriving Show
makeDomain :: Point
-> Double
-> Double
-> Double
-> 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
type DomainSeeds = (Seed, Seed, Seed, Seed, Seed, Seed)
getDimensions :: Domain -> (Double, Double, Double)
getDimensions (Domain xmin xmax ymin ymax zmin zmax) =
(xmax xmin, ymax ymin, zmax zmin)
getCenter :: Domain -> Point
getCenter (Domain xmin xmax ymin ymax zmin zmax) =
(xmin + (xmax xmin) / 2, ymin + (ymax ymin) / 2, zmin + (zmax zmin) / 2)
volume :: Domain -> Double
volume !(Domain xmin xmax ymin ymax zmin zmax) =
(xmax xmin) * (ymax ymin) * (zmax zmin)
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))
pureSpawnParticles :: Domain
-> Flow
-> Seed
-> (VU.Vector Particle, Seed)
pureSpawnParticles d flow s = purifyRandomST (spawnParticles d flow) s
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')
openBoundaryInjection :: DomainSeeds
-> Domain
-> Double
-> 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'))
clipToDomain :: Domain -> Ensemble -> Ensemble
clipToDomain (Domain xmin xmax ymin ymax zmin zmax) ens =
let
pred' :: Particle -> Bool
pred' !((x, y, z), _) =
xmax >= x && x >= xmin &&
ymax >= y && y >= ymin &&
zmax >= z && z >= zmin
in
filterEnsemble pred' ens
freeVolume :: Domain
-> Body
-> Int
-> 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)
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)