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
simulate :: Domain
-> Body
-> Flow
-> Time
-> Bool
-> Double
-> Double
-> Int
-> Surface
-> (Double, Double, Double)
-> Int
-> Int
-> IO (Int, Ensemble, MacroField)
simulate domain body flow
dt emptyStart ex sepsilon ssteps
surface
(mx, my, mz) volumePoints gsplit =
let
evolve :: (Ensemble, ParallelSeeds, DomainSeeds)
-> (Ensemble, ParallelSeeds, DomainSeeds)
evolve (ens, gseeds, dseeds) =
let
(e, dseeds') = openBoundaryInjection dseeds domain ex flow ens
(e', gseeds') = motion gseeds body dt surface e
e'' = clipToDomain domain e'
in
(e'', gseeds', dseeds')
macroSubdiv :: Grid
macroSubdiv = UniformGrid domain mx my mz
stabilized :: Ensemble -> Ensemble -> Bool
stabilized ens prevEns =
(abs $
((fromIntegral $ ensembleSize ens) /
(fromIntegral $ ensembleSize prevEns) 1)) < sepsilon
sim1 :: (Ensemble, ParallelSeeds, DomainSeeds)
-> Bool
-> Int
-> 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
gs <- replicateM gsplit $ randomSeed
(s1:s2:s3:s4:s5:s6:_) <- replicateM 6 randomSeed
vs <- replicateM gsplit $ randomSeed
startEnsemble <- if emptyStart
then return emptyEnsemble
else do
is <- randomSeed
return $ fst $ initializeParticles domain flow body is
fst <$> runMacroSampling
(sim1
(startEnsemble, gs, (s1, s2, s3, s4, s5, s6))
False 0)
vs macroSubdiv body volumePoints ssteps