
Percolation Simulation

This is a simulation of aspects of mathematical percolation theory.

For other percolation simulations, see Category:Percolation Simulations, or:
1D finite grid simulation
Mean run density 
2D finite grid simulations
Site percolation | Bond percolation | Mean cluster density

Given an 



M
×
N


{\displaystyle M\times N}

 rectangular array of cells numbered 




c
e
l
l

[
0..
M
−
1
,
0..
N
−
1
]


{\displaystyle \mathrm {cell} [0..M-1,0..N-1]}

, assume 



M


{\displaystyle M}

 is horizontal and 



N


{\displaystyle N}

 is downwards. Each 




c
e
l
l

[
m
,
n
]


{\displaystyle \mathrm {cell} [m,n]}

 is bounded by (horizontal) walls 




h
w
a
l
l

[
m
,
n
]


{\displaystyle \mathrm {hwall} [m,n]}

 and 




h
w
a
l
l

[
m
+
1
,
n
]


{\displaystyle \mathrm {hwall} [m+1,n]}

; (vertical) walls 




v
w
a
l
l

[
m
,
n
]


{\displaystyle \mathrm {vwall} [m,n]}

 and 




v
w
a
l
l

[
m
,
n
+
1
]


{\displaystyle \mathrm {vwall} [m,n+1]}



Assume that the probability of any wall being present is a constant 



p


{\displaystyle p}

 where

 



0.0
≤
p
≤
1.0


{\displaystyle 0.0\leq p\leq 1.0}


Except for the outer horizontal walls at 



m
=
0


{\displaystyle m=0}

 and 



m
=
M


{\displaystyle m=M}

 which are always present.

The task

Simulate pouring a fluid onto the top surface (



n
=
0


{\displaystyle n=0}

) where the fluid will enter any empty cell it is adjacent to if there is no wall between where it currently is and the cell on the other side of the (missing) wall.

The fluid does not move beyond the horizontal constraints of the grid.

The fluid may move “up” within the confines of the grid of cells. If the fluid reaches a bottom cell that has a missing bottom wall then the fluid can be said to 'drip' out the bottom at that point.

Given 



p


{\displaystyle p}

 repeat the percolation 



t


{\displaystyle t}

 times to estimate the proportion of times that the fluid can percolate to the bottom for any given 



p


{\displaystyle p}

.

Show how the probability of percolating through the random grid changes with 



p


{\displaystyle p}

 going from 



0.0


{\displaystyle 0.0}

 to 



1.0


{\displaystyle 1.0}

 in 



0.1


{\displaystyle 0.1}

 increments and with the number of repetitions to estimate the fraction at any given 



p


{\displaystyle p}

 as 



t
=
100


{\displaystyle t=100}

.

Use an 



M
=
10
,
N
=
10


{\displaystyle M=10,N=10}

 grid of cells for all cases.

Optionally depict fluid successfully percolating through a grid graphically.

Show all output on this page.




{-# LANGUAGE OverloadedStrings #-}
import Control.Monad
import Control.Monad.Random
import Data.Array.Unboxed
import Data.List
import Formatting
 
data Field = Field { f :: UArray (Int, Int) Char  
                   , hWall :: UArray (Int, Int) Bool  
                   , vWall :: UArray (Int, Int) Bool  
                   }
 
-- Start percolating some seepage through a field.
-- Recurse to continue percolation with new seepage.
percolateR :: [(Int, Int)] -> Field -> (Field, [(Int,Int)])
percolateR [] (Field f h v) = (Field f h v, [])
percolateR seep (Field f h v) = 
    let ((xLo,yLo),(xHi,yHi)) = bounds f
        validSeep = filter (\p@(x,y) ->    x >= xLo 
                                        && x <= xHi 
                                        && y >= yLo 
                                        && y <= yHi 
                                        && f!p == ' ') $ nub $ sort seep
 
        north (x,y) = if v ! (x  ,y  ) then [] else [(x  ,y-1)]
        south (x,y) = if v ! (x  ,y+1) then [] else [(x  ,y+1)]
        west  (x,y) = if h ! (x  ,y  ) then [] else [(x-1,y  )]
        east  (x,y) = if h ! (x+1,y  ) then [] else [(x+1,y  )]
        neighbors (x,y) = north(x,y) ++ south(x,y) ++ west(x,y) ++ east(x,y)
 
    in  percolateR 
            (concatMap neighbors validSeep)
            (Field (f // map (\p -> (p,'.')) validSeep) h v) 
 
-- Percolate a field;  Return the percolated field.
percolate :: Field -> Field
percolate start@(Field f _ _) = 
    let ((_,_),(xHi,_)) = bounds f
        (final, _) = percolateR [(x,0) | x <- [0..xHi]] start
    in final
 
-- Generate a random field.
initField :: Int -> Int -> Double -> Rand StdGen Field
initField width height threshold = do
    let f = listArray ((0,0), (width-1, height-1)) $ repeat ' '
 
    hrnd <- fmap (<threshold) <$> getRandoms
    let h0 = listArray ((0,0),(width, height-1)) hrnd       
        h1 = h0 // [((0,y), True) | y <- [0..height-1]]     -- close left 
        h2 = h1 // [((width,y), True) | y <- [0..height-1]] -- close right 
 
    vrnd <- fmap (<threshold) <$> getRandoms 
    let v0 = listArray ((0,0),(width-1, height)) vrnd   
        v1 = v0 // [((x,0), True) | x <- [0..width-1]]  -- close top
 
    return $ Field f h2 v1
 
-- Assess whether or not percolation reached bottom of field.
leaks :: Field -> [Bool]
leaks (Field f _ v) = 
    let ((xLo,_),(xHi,yHi)) = bounds f
    in [f!(x,yHi)=='.' && not (v!(x,yHi+1)) | x <- [xLo..xHi]]
 
-- Run test once; Return bool indicating success or failure.
oneTest :: Int -> Int -> Double -> Rand StdGen Bool
oneTest width height threshold = 
    or.leaks.percolate <$> initField width height threshold
 
-- Run test multple times; Return the number of tests that pass.
multiTest :: Int -> Int -> Int -> Double -> Rand StdGen Double
multiTest testCount width height threshold = do
    results <- replicateM testCount $ oneTest width height threshold
    let leakyCount = length $ filter id results
    return $ fromIntegral leakyCount / fromIntegral testCount
 
-- Helper function for display
alternate :: [a] -> [a] -> [a]
alternate [] _ = []
alternate (a:as) bs = a : alternate bs as
 
-- Display a field with walls and leaks.
showField :: Field -> IO ()
showField field@(Field a h v) =  do
    let ((xLo,yLo),(xHi,yHi)) = bounds a
        fLines =  [ [ a!(x,y) | x <- [xLo..xHi]] | y <- [yLo..yHi]]
        hLines =  [ [ if h!(x,y) then '|' else ' ' | x <- [xLo..xHi+1]] | y <- [yLo..yHi]]
        vLines =  [ [ if v!(x,y) then '-' else ' ' | x <- [xLo..xHi]] | y <- [yLo..yHi+1]]
        lattice =  [ [ '+' | x <- [xLo..xHi+1]] | y <- [yLo..yHi+1]]
 
        hDrawn = zipWith alternate hLines fLines
        vDrawn = zipWith alternate lattice vLines
    mapM_ putStrLn $ alternate vDrawn hDrawn
 
    let leakLine = [ if l then '.' else ' ' | l <- leaks field]
    putStrLn $ alternate (repeat ' ') leakLine
 
main :: IO ()
main = do
  g <- getStdGen
  let threshold = 0.45
      (startField, g2) = runRand (initField 10 10 threshold) g
 
  putStrLn ("Unpercolated field with " ++ show threshold ++ " threshold.")
  putStrLn ""
  showField startField
 
  putStrLn ""
  putStrLn "Same field after percolation."
  putStrLn ""
  showField $ percolate startField
 
  let testCount = 10000
      densityCount = 10
  putStrLn ""
  putStrLn ("Results of running percolation test " ++ show testCount ++ " times with thresholds ranging from 0/" ++ show densityCount ++ " to " ++ show densityCount ++ "/" ++ show densityCount ++ " .")
  let densities = [0..densityCount]
  let tests = sequence [multiTest testCount 10 10 v 
                           | density <- densities,
                             let v = fromIntegral density / fromIntegral densityCount ]
  let results = zip densities (evalRand tests g2)
  mapM_ print [format ("p=" % int % "/" % int % " -> " % fixed 4) density densityCount x | (density,x) <- results]