Simulating Experiences: Excursions in Programmming

Interfacial Diffusion Fronts and Gradient Percolation

by

Richard J. Gaylord and Kazume Nishidate

Published in
Mathematica in Education and Research
Vol.4 No.1
Copyright 1995 TELOS/Springer-Verlag

Introduction

We will show here how to determine the diffusion front. From this, various properties of the diffusion front, such as its fractal dimension, and its mean position and width can be computed.

The Gradient Percolation System

We will consider the diffusion of particles A on a two-dimensional lattice B with the diffusion source being a line of A atoms kept a constant concentration of one. The concentration of A particles will vary as a function of time and the distance from the source.

The diffusion front for the diffusion of A particles onto a lattice B over time can be generated in a straightforward manner (the interacting random walkers cellular automaton [Gaylord and Nishidate 1995] can be used to do this) as follows: a square lattice which is initially empty (all of the sites have value 0), evolves over many time steps in each of which the first row of the lattice is filled with particles (the values of all the sites in the first row are set to 1) and then all of the particles on the lattice try to move to empty nearest neighbor sites inside the lattice.

Rather than generating a diffusion front in this computationally intensive way, we can simply create a lattice with a gradient of concentration of occupied sites. This is done by distributing particles randomly on each row of the lattice with a density distribution probability, p, that decreases monotonically with the distance of a row from the first row. The precise form of p is not important (eg., the concentration of particles can have a linear or erfc dependence on position) as long as it fills the first row and leaves the last row empty.

Note that while distributing occupied sites on a lattice in this way is characteristic of percolation [see chpt. 5 in Gaylord and Wellin 1995], the clusters that are generated in this model differ from the clusters generated in the usual percolation system in that the probabiliy of a site being occupied in ordinary percolation is constant over the whole lattice (ie., the same for all sites in the lattice). The directional dependence in the system being considered here is accordingly referred to as gradient percolation.

The system consists of a rectangular lattice having m rows and n columns.  The probability of a site being occupied in a given row decreases linearly with the distance of the row from the first row

The lattice can be generated using the following program:

gradientPercolation[n_, m_] :=
Module[{p = 0},
   createRow = 
   (p += 1/(n - 1); h[#, Table[Floor[Random[] + p], {m}]])&; 
  Apply[List, Flatten[Nest[createRow, Table[0, {m}], n - 1],
        Infinity, h] ]
       ]

A typical gradient percolation system is shown below.

Show[Graphics[
   RasterArray[gradientPercolation[100, 100] /. 
       {0 -> GrayLevel[0], 1 -> GrayLevel[1]}]
             ],
       AspectRatio->Automatic]

Identifying the Diffusion Front

A diffusion front (or 'seashore' or 'frontier') is defined by looking at the cluster of occupied (or A) sites containing the occupied sites in last row of the lattice and the cluster of empty (or B) sites containing the empty sites in the first row. The interface consisting of the sites which  give the closest possible contact between these two 'infinite' clusters can be obtained by determining the sites in these two clusters in different ways:

The sites in the cluster of occupied sites are in the von Neumann neighborhood of other sites in the cluster (the von Neumann neighborhood consists of the four nearest neighbor sites lying N, E, S, W of a site).

The sites in the cluster of empty sites are in the Moore neighborhood of other sites in the cluster  (the Moore neighborhood consists of the eight nearest neighbor sites lying N, NE, E, SE, S, SW, W, NW of a site) .

Two diffusion fronts can be defined:

(1) the sites in the 'infinite' empty site cluster that lie in the von Neumann neighborhood of a site in the 'infinite' occupied site cluster.

(2) the sites in the 'infinite' occupied site cluster that lie in the Moore neighborhood of a site in the 'infinite' empty site cluster.

In determining the diffusion fronts, the lattice is taken to have periodic boundary condtions in the direction perpendicular to the concentration gradient.

We will show how to determine the diffusion front of a two-dimensional gradient percolation system in a step-by-step fashion. Below is the simple system which we will use to illustrate the computation. The occupied sites have value 1 and the unoccupied sites have value 0.

sapoval = 
 {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
  {0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0}, 
  {0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0}, 
  {0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0}, 
  {1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0}, 
  {0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0}, 
  {0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0}, 
  {1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0}, 
  {1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, 1, 0, 1}, 
  {1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1}, 
  {1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0}, 
  {1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0}, 
  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1}, 
  {1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1}, 
  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}, 
  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}};

We first identify the clusters of occupied sites in the sapoval lattice. This can be done using the following cluster labeling program [see p. 154 in Gaylord and Wellin 1995].

clusterLabelvonNeumann[lat_List]:=
  Module[{i=2, clusterCornerID, cornerLabels, reLabel},
        clusterCornerID[1,0, 0]:= i++;
        clusterCornerID[a_, __]:= a;
        Attributes[clusterCornerID] = Listable;

        cornerLabels = 
          clusterCornerID[#, 
                          RotateRight[#,{1,0}],
                          RotateRight[#,{0,1}]]&[lat];

        reLabel[0, ___] := 0;
        reLabel[a_, b_, c_, d_, e_] := Max[a, b, c, d, e];
        Attributes[reLabel] = Listable;

        FixedPoint[reLabel[#,
                           RotateRight[#,{1, 0}],
                           RotateRight[#,{0, 1}],
                           RotateRight[#,{-1,0}],
                           RotateRight[#,{0,-1}]]&,
                   cornerLabels]
        ]

Applying the clusterLabelvonNeumann function to the sapoval lattice and then re-numbering the identified clusters so that their numbers run sequentially with no gaps, we get:

AClustersLat = clusterLabelvonNeumann[sapoval]; 

AClusters =
 AClustersLat /. 
  MapThread[Rule,
            {Reverse[Rest[Union[Flatten[AClustersLat]]]],
             Range[Length[Union[Flatten[AClustersLat]]] - 1]}
          ]

The clusterLabelvonNeumann program labels clusters of connected occupied sites (ie., sites that have value 1 and lie in the vonNeumann neighborhood of other sites that have value 1), starting in the bottom row and proceeding upwards. Therefore the sites having value 1 in AClusters belong to the 'infinite' occupied site cluster. Their locations are determined using the Position function.

infiniteAClusterlocs = Position[AClusters, 1]

To identify the clusters of unoccupied sites, we modify the clusterLabelvonNeumann program so that it labels connected sites that have value 1 and lie in the Moore neighborhood of other sites that have value 1.

clusterLabelMoore[lat_List]:=
  Module[{i = 2, clusterCornerID, cornerLabels, reLabel},
        clusterCornerID[1, 0, 0, 0] := i++;
        clusterCornerID[a_, __] := a;
        Attributes[clusterCornerID] = Listable;

        cornerLabels = 
          clusterCornerID[#, 
                          RotateRight[#, {1, 0}],
                          RotateRight[#, {0, 1}],
                          RotateRight[#, {1, 1}]]&[lat];

        reLabel[0, ___] := 0;
        reLabel[a_, b_, c_, d_, e_, f_, g_, h_, i_] := 
                          Max[a, b, c, d, e, f, g, h, i];
        Attributes[reLabel] = Listable;
        
        FixedPoint[
           reLabel[#,
                    RotateRight[#,{1, 0}],
                    RotateRight[#,{0, 1}],
                    RotateRight[#,{-1,0}],
                    RotateRight[#,{0,-1}],
                    RotateRight[#,{1, 1}],
                    RotateRight[#,{1, -1}],
                    RotateRight[#,{-1, 1}],
                    RotateRight[#,{-1,-1}]]&, 
           cornerLabels]
        ]

The clusterLabelMoore program labels clusters in the same way as the clusterLabelvonNeumann program, starting in the bottom row. Therefore, in order to use the program to identify the empty site clusters in the sapoval lattice, we change the 0's into 1's and 1's into 0's and reverse the order of the rows in the lattice.

BClustersLat = 
  clusterLabelMoore[Reverse[sapoval /. {0 -> 1, 1 -> 0}]]

Taking BClustersLat, reversing the rows and renumbering the clusters of empty sites, we get

BClusters = 
 Reverse[BClustersLat /. 
  MapThread[Rule,
           {Reverse[Rest[Union[Flatten[BClustersLat]]]],
            Range[Length[Union[Flatten[BClustersLat]]] - 1]}]
           ]

The sites having value 1 in BClusters belong to the infinite empty site cluster in the sapoval lattice. Their locations in the lattice is given by

infiniteBClusterlocs = Position[BClusters, 1]

Now that we have the list, infiniteAClusterlocs, of the sites in the infinite occupied site cluster, and the list, infiniteBClusterlocs, of the sites in the infinite empty site cluster, it is straightforward to determine the list of sites, frontierBSites, in infiniteBClusterlocs that lie in the von Neumann neighborhood of sites in infiniteAClusterlocs, and the list of sites, frontierASites, in infiniteAClusterlocs that lie in the Moore neighborhood of sites in infiniteBClusterlocs.

frontierBSites = 
Select[infiniteBClusterlocs, 
MemberQ[
  Flatten[
   Map[Function[y,Map[(y + #)&,{{1,0},{-1,0},{0,1},{0,-1}}]],
    infiniteAClusterlocs], 1
         ], #
       ]&]

frontierASites = 
Select[infiniteAClusterlocs, 
MemberQ[
 Flatten[
 Map[Function[y,
   Map[(y + #)&,
     {{1,0},{-1,0},{0,1},{0,-1},{1,1},{1,-1},{-1,1},{-1,-1}}]
              ],
      infiniteBClusterlocs
    ], 1], #]&
     ]

Using the sapoval lattice and the lists of the sites along the two diffusion fronts, frontierASites and frontierBSites, we can create the a graphic in which the A sites are colored blue, the  A frontier sites are colored green, and the B frontier sites are colored red.

Show[Graphics[RasterArray[
  Reverse[ReplacePart[ReplacePart[sapoval, FA, frontierASites], FB, frontierBSites]] /. 
     {0 -> GrayLevel[1], 1 -> GrayLevel[0], 
      FA -> GrayLevel[.33], FB -> GrayLevel[.67]}
                        ], 
  AspectRatio -> Automatic, Frame -> True, 
  FrameTicks -> None]
    ]

Computer Simulation Projects

1. Create a program for calculating the diffusion fronts of a two-dimensional gradient percolation system, using the code fragments presented here.

2. An analytical  analysis indicates that the concentration, p[x] follows  so-called Fickian behavior, going as Erfc[x/L] where x is the distance from the source and L = 2 (D*t)^0.5 is the diffusion length at time t with diffusion coefficient D. Using this concentration profile, determine the power law dependences of the number of particles on the front and the spread (width) of the front on the diffusion length.

3. A three-dimensional gradient percolation system can be created using

gradientPercolation3D[n_, m_] :=
Module[{p = 0}, 
  laydown =
     (p += 1/(n - 1); 
      Append[#, Table[Floor[Random[] + p], {m}, {m}]]
     )&;

   Nest[laydown,{Table[0, {m}, {m}]}, n - 1]
     ]

A graphic of the occupied sites in the system can be created by locating the positions of sites with value 1 in the lattice and applying the Cuboid graphic primitive to those positions. 

For example,

Show[Graphics3D[Map[Cuboid, 
      Position[gradientPercolation3D[9, 10], 1]]]]

Defining the interface so that the boundaries of the two 'infinite' clusters are in contact with each other everywhere, each site in the cluster of occupied sites is one of the four nearest neighbor sites in the  N, E, S, W directions of another site in the cluster, while each site in the cluster of empty sites is one of the 12 sites that are either nearest neighbor sites in the N, NE, E, SE, S, SW, W, NW directions, or next nearest neighbor sites in the N, E, S, W directions of another site in the cluster.

General References

Gaylord, Richard J. and Wellin, Paul R. "Computer Simulations with Mathematica: Explorations in Complex Physical and Biological Systems", TELOS/Springer-Verlag (1995) Chapter ABC

Sapoval, B., Rosso, M., and Gouyet, J. F. "The fractal nature of a diffusion front and the relation to percolation". J. Physique Lett., 46 (1985), L149-L156

Rosso, M., Gouyet, J. F., and Sapoval, B. "Gradient Percolation in Three Dimensions and Relation to Diffusion Fronts". Phys. Rev. Lett., 57 (1986), 3195-3198.

Wool, R. P. and Long, J. M.. "Fractal Nature of Polymer Interfaces". Macromolecules, 26 (1993), 5227-5239.