Permalink
Browse files

Mssive update to hg default 756.

  • Loading branch information...
1 parent cbabbd8 commit a49f1c0f0917b6debf386436eb34521cee19df27 Timothy Kelley committed Mar 30, 2012
View
@@ -0,0 +1,14 @@
+import qualified System.Random.MWC as MWC
+import System.Random.LFG (defaultLags, largeLag, smallLag)
+import System.Random.LFG.Pure as LFG
+import Control.Monad
+
+main =
+ do
+ mwc <- MWC.create
+ initials <- replicateM (3 * largeLag defaultLags) (MWC.uniform mwc)
+ print $ sum initials
+ print $ smallLag defaultLags
+ let rngs = LFG.sequences defaultLags initials
+ print $ (rngs !! 2) !! 417
+
@@ -0,0 +1,75 @@
+-- Histogram_Test.hs
+-- T. M. Kelley
+-- Sep 02, 2011
+-- (c) Copyright 2011 LANSLLC, all rights reserved
+
+module Test.Histogram_Test (tests) where
+
+import Test.Framework (testGroup)
+import Test.Framework.Providers.HUnit
+import Test.Framework.Providers.QuickCheck2 (testProperty)
+import Test.HUnit
+import Test.QuickCheck
+import Test.Arbitraries
+
+import Histogram
+import Physical
+import SoftEquiv
+
+import qualified Data.Vector.Unboxed as U
+
+-- * count
+
+-- | count with a null weight should leave the counts & squares unchanged
+prop_nullCountUnchanged :: HistNEnergy -> Bool
+prop_nullCountUnchanged (HistNEnergy h e) = cmpEHists h' h (EnergyWeight 1e-16)
+ where h' = count h (e,EnergyWeight 0)
+
+-- | count is commutative to within FP error (on same index)
+-- (that is, {count h e1; count h e2} == {count h e2; count h e1})
+prop_cmtvSameIdx :: HistNEnergy -> EnergyWeight -> EnergyWeight -> Bool
+prop_cmtvSameIdx (HistNEnergy h e) ew1 ew2 = cmpEHists h1 h2 (EnergyWeight 1e-14)
+ where h1 = count (count h (e,ew1)) (e,ew2)
+ h2 = count (count h (e,ew2)) (e,ew1)
+
+cmpEHists ha hb tol = U.foldl1' (&&) compCounts
+ where compCounts = U.zipWith cmpfunc (ehcounts ha) (ehcounts hb)
+ compSquares = U.zipWith cmpfunc (ehsquares ha) (ehsquares hb)
+ cmpfunc x y = softEquiv x y tol
+
+-- | count is absolutely commutative (on different indices)
+
+
+-- findBin properties
+-- bb ! (findBin e bb) <= e
+prop_fbLowerBound (HistNEnergy (EHist _ _ _ bbs) e) = bbs U.! i <= e
+ where i = findBin e bbs
+
+-- bb ! (findBin e bb) + 1 > e
+prop_fbUpperBound (HistNEnergy (EHist _ _ _ bbs) e) = bbs U.! (i + 1) > e
+ where i = findBin e bbs
+
+-- throws exception if asked for an energy outside of the bin boundaries
+
+-- divIfne0
+prop_divIfne0_finite :: Double -> Double -> Bool
+prop_divIfne0_finite x y = q /= 1/0 && q /= -1/0
+ where q = x `divIfne0` y
+
+
+-- aggregate tests for framework
+tests = [testGroup "histogram tests"
+ [
+ testProperty "count 0 wt: histogram unchanged" prop_nullCountUnchanged
+ ,testProperty "count commutes (within FP on same idx)" prop_cmtvSameIdx
+ -- ,testProperty "count commutes (absolutely on different idx)" prop_cmtvDiffIdx
+ ,testProperty "findBin is lower bound on e" prop_fbLowerBound
+ ,testProperty "findBin + 1 is upper bound on e" prop_fbLowerBound
+ ]
+ ]
+
+
+-- version
+-- $Id$
+
+-- End of file
@@ -0,0 +1,226 @@
+-- PickEvent_Test.hs
+-- T. M. Kelley
+-- Aug 08, 2011
+-- (c) Copyright 2011 LANSLLC, all rights reserved
+
+module Test.PickEvent_Test (tests) where
+
+import Test.Framework (testGroup)
+import Test.Framework.Providers.HUnit
+import Test.Framework.Providers.QuickCheck2 (testProperty)
+import Test.HUnit
+import Test.QuickCheck
+import Test.Arbitraries
+
+-- The library under test
+import Sphere1D
+import Mesh
+import Cell
+import Physical
+import PRNG
+import Material
+import SoftEquiv
+import MC
+import Sigma_HBFC
+import Particle as P
+import Constants
+
+-- helpers
+import qualified Data.Vector as V
+import Data.List (zip4,unfoldr)
+
+import Control.Applicative ()
+
+-- | exports from this module
+tests = [ testGroup "pickEvent" pECases]
+
+
+-- * pickEvent
+pECases = [ pECase1
+ , pECase2
+ , pECase3
+ , pECase4
+ , pECase5
+ ]
+pECase1 = testCase "pickEvent: stream 0.5 cm to high boundary" pETest1
+pECase2 = testCase "pickEvent: stream through 10 cells, escape on last one"
+ pETest2
+pECase3 = testCase "pickEvent: Nucleon absorption" pETest3
+pECase4 = testCase "pickEvent: Nucleon elastic scatter" pETest4
+pECase5 = testCase "pickEvent: e- scatter" pETest5
+
+-- These test cases changed quite a bit upon refactoring
+-- MC. I've left a bunch of the old expressions here--they should
+-- be incorporated into new cases that test the new parts of MC.
+
+pETest1 :: Assertion
+pETest1 =
+ let
+ msh = genMesh 10 1.0
+ rng = makeFakeGen $ cycle [0.30897681610609407,
+ 0.92436920169905545,
+ 0.21932404923057958]
+ p = Particle (Position 0.5) (Direction 1) (Time 1) (Energy 5) (EnergyWeight 1) (CellIdx 0)
+ event = runRnd rng (pickEvent nuE msh p)
+ eventExp = BoundaryCand (Distance 0.5) Hi
+ -- eventExp = Boundary Transmit (Distance 0.5) Hi (Energy 5) (EnergyWeight 1)
+ in assertBool ("chosen event, " ++ show event ++
+ ", did not match expected event, " ++
+ show eventExp)
+ (event == eventExp)
+
+{- This case runs a particle through a 10 cell mesh. On the last step, -
+ - it should escape. Not the most convincing test in the world, meh. -}
+pETest2 :: Assertion
+pETest2 =
+ let
+ msh = genMesh 10 1.0
+ rng = makeFakeGen $ cycle [0]
+ p0 = Particle (Position 0.5) (Direction 1) (Time 1) (Energy 5) (EnergyWeight 1) (CellIdx 0)
+ ps = take 10 (unfoldr pf p0)
+ -- pf: evolve particle state
+ pf :: Particle -> Maybe (Particle,Particle)
+ pf p@(Particle {P.pos = Position x,cellIdx = CellIdx cidx}) = Just $
+ (p, p{ P.pos = Position (x + d), cellIdx = CellIdx (cidx + 1)})
+ where (Distance d) = candDist evt
+ evt = runRnd rng (pickEvent nuE msh p)
+ p10 = last ps
+ event = runRnd rng (pickEvent nuE msh p10)
+ eventExp = BoundaryCand (Distance 1.0) Hi
+ -- eventExp = Boundary Escape (Distance 1.0) Hi (Energy 5) (EnergyWeight 1)
+ in assertBool ("chosen event, " ++ show event ++
+ ", did not match expected event, " ++
+ show eventExp)
+ (event == eventExp)
+
+{- runs a particle through a cell with reasonable nucleon density; RNG -
+ - stacked to select nucleon absorption. Note that the sequence of RN -
+ - consumption is different from the Python/C++ version. -}
+pETest3 :: Assertion
+pETest3 =
+ let
+ n = 10
+ dx = 1e6
+ clls = genCells (genBndConds n) (genBounds n dx) denseNMat
+ msh = Sphere1D . V.fromList $ clls
+ rng = makeFakeGen $ cycle [0.30897681610609407, -- this is used for d_coll
+ -- 0.92436920169905545, -- not used in Haskell
+ 0.21932404923057958]
+ p = Particle (Position 0.5) (Direction 1) (Time 100) (Energy 5) (EnergyWeight 1) (CellIdx 0)
+ event = runRnd rng (pickEvent nuE msh p)
+ dExp = 7343.827
+ -- mExp = 5/c
+ -- eExp = 5
+ eventExp = CollisionCand (Distance dExp)
+ d = Physical.distance . candDist $ event
+ -- eventExp = Collision NuclAbs (Distance dExp) (Momentum mExp) (Energy eExp)
+ in assertBool ("chosen event, " ++ show event ++
+ ", did not match expected event, " ++
+ show eventExp)
+ (softEquiv d dExp 1e-7)
+ -- (case event of
+ -- Collision NuclAbs (Distance d) (Momentum pdep) (Energy edep) ->
+ -- softEquiv d dExp 1e-3 &&
+ -- softEquiv pdep mExp 1e-10 &&
+ -- eExp == edep
+ -- _ -> False
+ -- )
+
+{- runs a particle through a cell with reasonable nucleon density; RNG -
+ - stacked to select nucleon scattering. Note that the sequence of RN -
+ - consumption is different from the Python/C++ version. -}
+pETest4 :: Assertion
+pETest4 =
+ let
+ n = 10
+ dx = 1e6
+ clls = genCells (genBndConds n) (genBounds n dx) denseNMat
+ msh = Sphere1D . V.fromList $ clls
+ rng = makeFakeGen $ cycle [0.21932404923057958, -- used for d_coll
+ -- 0.20867489035315723,
+ 0.91525579001682567]
+ p = Particle (Position 0.5) (Direction 1) (Time 100) (Energy 5) (EnergyWeight 1) (CellIdx 0)
+ event = runRnd rng (pickEvent nuE msh p)
+ dExp = 9486.756524
+ -- mExp = 5/c
+ -- eExp = 0
+ eventExp = CollisionCand (Distance dExp)
+ -- eventExp = Collision NuclEl (Distance dExp) (Momentum mExp) (Energy eExp)
+ d = Physical.distance . candDist $ event
+ in assertBool ("chosen event, " ++ show event ++
+ ", did not match expected event, " ++
+ show eventExp)
+ (softEquiv d dExp 1e-9)
+ -- (case event of
+ -- Collision NuclEl (Distance d) (Momentum pdep) (Energy edep) ->
+ -- softEquiv d dExp 1e-3 &&
+ -- -- softEquiv pdep mExp 1e-10 &&
+ -- eExp == edep
+ -- _ -> False
+ -- )
+
+{- runs a particle through a cell with reasonable lepton density; RNG -
+ - stacked to select lepton scattering. Note that the sequence of RN -
+ - consumption is different from the Python/C++ version. -}
+pETest5 :: Assertion
+pETest5 =
+ let
+ n = 10
+ dx = 1e6
+ clls = genCells (genBndConds n) (genBounds n dx) denseEMinusMat
+ msh = Sphere1D . V.fromList $ clls
+ rng = makeFakeGen $ cycle [0.9, -- used for d_coll
+ -- 0.9,
+ 0.1]
+ p = Particle (Position 0.5) (Direction 1) (Time 100) (Energy 5) (EnergyWeight 1) (CellIdx 0)
+ event = runRnd rng (pickEvent nuE msh p)
+ dExp = 38310.45776
+ -- mExp = 1/c
+ -- eExp = 1
+ eventExp = CollisionCand (Distance dExp)
+ -- eventExp = Collision EMinusInel (Distance dExp) (Momentum mExp) (Energy eExp)
+ d = Physical.distance . candDist $ event
+ in assertBool ("chosen event, " ++ show event ++
+ ", did not match expected event, " ++
+ show eventExp)
+ (softEquiv d dExp 1e-9)
+ -- (case event of
+ -- Collision EMinusInel (Distance d) (Momentum pdep) (Energy edep) ->
+ -- softEquiv d dExp 1e-3 &&
+ -- -- softEquiv pdep mExp 1e-10 &&
+ -- eExp == edep
+ -- _ -> False
+ -- )
+
+-- helpers
+genBndConds :: Int -> [BoundaryCondition]
+genBndConds n = [Refl] ++ [Transp | i<-[1..n-1]] ++ [Vac]
+
+genBounds :: Int -> Double -> [Position]
+genBounds n dx = [Position $ dx * fromIntegral i | i <- [0..n]]
+
+zeroMat :: Material
+zeroMat = Material (Velocity 0) (Temperature 0)
+ (Density 1.0) (NDensity 0) (NDensity 0)
+
+denseNMat :: Material
+denseNMat = Material (Velocity 0) (Temperature 0)
+ (Density 1e14) (NDensity 0) (NDensity 0)
+
+denseEMinusMat :: Material
+denseEMinusMat = Material (Velocity 0) (Temperature 1)
+ (Density 0) (NDensity (1e14/pmg)) (NDensity 0)
+
+genCells :: [BoundaryCondition] -> [Position] -> Material -> [Cell]
+genCells bcs xs matl = [Cell posl posh bcl bch matl |
+ (posl,posh,bcl,bch) <- zip4 (init xs)
+ (tail xs)
+ (init bcs)
+ (tail bcs)]
+
+genMesh :: Int -> Double -> Sphere1D
+genMesh n dx = Sphere1D . V.fromList $
+ genCells (genBndConds n) (genBounds n dx) zeroMat
+
+
+-- end of file
Oops, something went wrong.

0 comments on commit a49f1c0

Please sign in to comment.