
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.

Assume that the probability of any cell being filled is a constant 



p


{\displaystyle p}

 where

 



0.0
≤
p
≤
1.0


{\displaystyle 0.0\leq p\leq 1.0}


The task

Simulate creating the array of cells with probability 



p


{\displaystyle p}

 and then
testing if there is a route through adjacent filled cells from any on row 



0


{\displaystyle 0}

 to any on row 



N


{\displaystyle N}

, i.e. testing for site percolation.

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
=
15
,
N
=
15


{\displaystyle M=15,N=15}

 grid of cells for all cases.

Optionally depict a percolation through a cell grid graphically.

Show all output on this page.

import std.stdio, std.random, std.array, std.datetime;
 
enum size_t nCols = 15,
            nRows = 15,
            nSteps = 11,     // Probability granularity.
            nTries = 20_000; // Simulation tries.
 
enum Cell : char { empty   = ' ', filled  = '#', visited = '.' }
alias Grid = Cell[nCols][nRows];
 
void initialize(ref Grid grid, in double probability, ref Xorshift rng) {
    foreach (ref row; grid)
        foreach (ref cell; row)
            cell = (rng.uniform01 < probability) ? Cell.empty : Cell.filled;
}
 
void show(in ref Grid grid) @safe {
    writefln("%(|%(%c%)|\n%)|", grid);
}
 
bool percolate(ref Grid grid) pure nothrow @safe @nogc {
    bool walk(in size_t r, in size_t c) nothrow @safe @nogc {
        enum bottom = nRows - 1;
        grid[r][c] = Cell.visited;
 
        if (r < bottom && grid[r + 1][c] == Cell.empty) { // Down.
            if (walk(r + 1, c))
                return true;
        } else if (r == bottom)
            return true;
 
        if (c && grid[r][c - 1] == Cell.empty) // Left.
            if (walk(r, c - 1))
                return true;
 
        if (c < nCols - 1 && grid[r][c + 1] == Cell.empty) // Right.
            if (walk(r, c + 1))
                return true;
 
        if (r && grid[r - 1][c] == Cell.empty) // Up.
            if (walk(r - 1, c))
                return true;
 
        return false;
    }
 
    enum startR = 0;
    foreach (immutable c; 0 .. nCols)
        if (grid[startR][c] == Cell.empty)
            if (walk(startR, c))
                return true;
    return false;
}
 
void main() {
    static struct Counter {
        double prob;
        size_t count;
    }
 
    StopWatch sw;
    sw.start;
 
    enum probabilityStep = 1.0 / (nSteps - 1);
    Counter[nSteps] counters;
    foreach (immutable i, ref co; counters)
        co.prob = i * probabilityStep;
 
    Grid grid;
    bool sampleShown = false;
    auto rng = Xorshift(unpredictableSeed);
 
    foreach (ref co; counters) {
        foreach (immutable _; 0 .. nTries) {
            grid.initialize(co.prob, rng);
            if (grid.percolate) {
                co.count++;
                if (!sampleShown) {
                    writefln("Percolating sample (%dx%d, probability =%5.2f):",
                             nCols, nRows, co.prob);
                    grid.show;
                    sampleShown = true;
                }
            }
        }
    }
    sw.stop;
 
    writefln("\nFraction of %d tries that percolate through:", nTries);
    foreach (const co; counters)
        writefln("%1.3f %1.3f", co.prob, co.count / double(nTries));
 
    writefln("\nSimulations and grid printing performed" ~
             " in %3.2f seconds.", sw.peek.msecs / 1000.0);
}