Diffusion-limited aggregation
simulations of fractals

1998 Scientific Computing III Project

Paul Titze <ptitze@vislab.usyd.edu.au>
Braden Kohary <bkohary@vislab.usyd.edu.au>

Abstract

Various shapes such as a bunch of grapes on an overhead vine, dendrites, colloids appear to arise from an aggregation process that is limited by diffusion. To investigate this process, various simulations of diffusion-limited aggregations are accomplished starting from a 2D lattice containing a seed particle in the middle. A 3D version of the simulation is also investigated and various observations are made on the fractal shapes of these clusters and their fractal dimensions are evaluated.

Keywords: Diffusion-limited aggregation, fractal shapes, diffusion processes, particle cluster growth

Introduction

This project further investigates diffusion-limited aggregation processes described in Landau[1] which showed to have interesting areas of further study. A 2D lattice containing a seed particle in the middle is first used and itís fractal dimension is evaluated. A similar procedure was then used with 3D fractal clusters with a wide range of iterations and their fractal dimensions evaluated using boxes instead of squares. The simulation was then ported to a large computer with multiprocessors with the number of iterations increased to several million.

Theoretical Background

A 2D lattice containing a seed particle in the middle is first used and a circle a drawn around the particle and another particle is placed on the circumference of the circle at some random angle between 0 and 2(, the x and y positions of this new particle on the circle are obtained. This second particle is then released and executes a random walk restricted to vertical or horizontal jumps between lattice sites and this is determined by generating a uniform random number 0 < rxy < 1 and applying the rule.

if rxy < 0.5 , the motion is vertical

if rxy > 0.5 , the motion is horizontal

This is a type of Brownian motion that simulates the diffusion process and has a fractal nature. The model is made more realistic by allowing the length for each step vary according to a random Gaussian distribution. If at some point during its random walk the particle finds another particle within one lattice spacing, they stick together and the walk terminates. If the particle passes outside the circle from which it was released, it is lost forever. The process is repeated as often as desired and a typical cluster growth shown in Fig 1. was obtained by running the code provided by Landau in DLA.C (Appendix A). Because many particles ìget lostî, it is necessary to generate hundreds of thousands of particles to form a cluster of several hundred particles. Most of the simulations show that the number of actual particles generated if of the order of a hundred less then the number of iterations used.

The number of iterations in our code has been kept to below 400000 iterations to keep the simulation time reasonably low. The coordinates for the various points are written to a file called DLA.DAT. A higher number of iterations of the order of a few million was used for the

The fortran version of our code, Appendix D.

Objects which have a fractal structure, are those which have the same structure in a small region as there is in the entire figure, in other words, if the object were infinitely dense, any part of the object could be scaled up in size and will be similar to the whole. This property is called self-similarity. If we continue to build the fractal object and look at how its densities (mass/area) changes with size, we see that as we continue the construction process, the density of each new structure is æ that of the previous one. This is unusual.

For ordinary objects the density is an intensive quantity independant of the size of an object. For this strange object, there is a power-law dependence of the desnity on the size of the object

Hence if the object in question has a fractional dimension of df , it has a fractal structure.

Obtaining the data and plots

The original code used as a starting point for this project was obtained from Landau[1] which gives a 2D simulation of the DLA process and is written in C, Appendix A. The code gives fractal shapes typical of Figure 1. To visualize the data written to DLA.DAT, GnuPlot was used to visualize the data for preliminary plots and to confirm the code was working correctly. A sample of a plot for 100000 iterations is shown in Appendix H. This code was then modified with extra code to compute the fractal dimension of the fractal 2D cluster. The code implements the procedure Landau used for the calculation by using several squares growing larger from the seed particle and calculating the density by knowing the number of particles within each given square. The length of the squares were increased with a unit length of 5 for each step. The densities for each given square was then calculated and the results, square length, density are written to a file for plotting in Mathematica. The code to implement the 2D fractal analysis can be found in Appendix B.

The generated data written to DLA.DAT was then imported into Mathematica for plotting and obtaining a log-log plot for the results. The slope from the fit of the points in the graph was then obtained and the fractal dimension of the cluster was then evaluated.

The fractal dimension df was then used to confirm whether our generated 2D and 3D clusters were a of fractal nature or not.

Working from the original DLA.C code from Landau, a 3D version of the simulation was coded in C to obtain a better visualization of the fractal shape of the diffusion-limited aggregation process which can be found in Appendix C, DLA3D.C. This code essentially implements the original 2D procedures in 3 Dimensions by choosing 2 random numbers along the surface of a sphere. Running the simulation for 100, 200, 300 and 400000 iterations, 3D DLA clusters were obtained and can be found in the results section. A number of software packages were used to visualize the data. GnuPlot showed to be inadequate for this purpose and initially used AVS, a more advanced package. An example of a 3D plot made with AVS is shown in Appendix G for 100000 iterations. AVS provided an adequate environment for plotting our data points but due to the extreme slow performance of the package, we had to move to Mathematica and plot the 3D clusters from there. For a 200000 iteration simulation, which corresponded to 5695 points in the DLA fractal cluster, AVS was not capable of handling this amount of points while Mathematica could handle a 400000 iteration run which corresponded to 37949 points in DLA.DAT.

Continuing with the implementation of our 3D code, the fractal analysis component was added which allowed the computation of the fractal dimension of the 3D cluster using a similar approach to the 2D version, the code is shown in Appendix E.

DLA3D.C was ported to a Fortran90 version, DLA3D.F90 (Appendix D) with the intention of running these simulations for several million iterations and to obtain larger fractal shapes in 3D. The Fortran90 version of our code was then optimized to take advantage of the parallel processing capabilities of the 20 processor Sydney Vislab supercomputer, Napier. This allowed us to gain first hand experience using high performance computing technics for parallel processors and learned how to use High Performance Fortran, HPF necessary to run our code on Napier.

The simulations that were run on Napier using our Fortran90 code and HPF are shown in Appendix I. The fortran code on Napier was slightly altered to display the number of processors used during the run. The procedure for running the code on Napier is slightly different to directly running A.OUT on the Silicon Graphics machines used for most of our work. A script is created which invokes our code, sets various environment variables such as the number of processors we wish to use, starts HPF to compile our Fortran90 code and submits it a Napier queue for execution. Our simulation is given a job number similar to printing queues and then executed. The main job queue on Napier which allowed the highest memory usage had a number of jobs running so our simulations used the ìNormal Queueî. On the completition of our run, Napier Emailed us with statistics of the run which are also shown in Appendix I. Note the large number of particles generated in the aggregation. A number of runs were made for 5, 5, 50, 50 million iterations with 4, 14, 4, 20 processors respectively. The generated DLA.DAT was too large to plot in any of the packages we used. However, the whole exercise on Napier prooved to be very enlighting and we now have a good understanding using Napier for large computation projects.

While plotting our points to obtain 3D fractal clusters using Mathematica, we investigated if there was some software available which could show our models in 3D in real time such that we could turn the 3D cluster at will and gain a better view of our cluster from various sides and angles. Although AVS showed that have some capabilities in this regard, it was too slow to be of practical use (this was using an O2 Silicon Graphics machine). We turned our attention to using Virtual Reality Modeling Language. VRML is a language that allows complete interactive virtual worlds to be experienced in full 3D with 3D spatial sound and animation through the Web. Using the generated data from DLA3D.DAT, VRML models of our fractal clusters were produced and viewed through Netscape. The sample code of DLA3D.WRL can be found in Appendix F.

Results

Using DLA3D.DAT, the fractal clusters were plotted in Mathematica by reading in the coordinate points for each particle then doing a 3Dplot. Several clusters were obtained for 100 to 400000 iteration simulations for both the C and Fortran90 versions of our code. Figure 2 shows a DLA cluster obtained for 100000 iterations which resulted in 1119 particles generated and 758 particles for the Fortran version.

Figure 2 - 100K iterations for dla3d.c

The fortran90 version of our code was also used to obtain the 3D clusters and gave a check on whether the code was correctly implementing the C version. Figure 3 shows a 100K run for the fortran90 version.

Figure 3 - 100K iterations for dla3d.f90

Several more runs were made with a higher number of iterations. Figures 4 and 5 show the resulting clusters obtained for 200K iterations for the C and Fortran90 versions respectively. Note the centre region towards the seed particle which is starting to become very dense. Throughout all these runs, the seed particle is located at (50,50,50).

Figure 4 - 200K iterations for dla3d.c

The fortran version shown in figure 5 resulted in 4452 particles generated. The process was repeated again for 300K iterations and figures 6 and 7 show the resulting clusters for the C and Fortran version respectively. The C version resulted in 15383 particles and the fortran version generated 14686 particles. Distinguishing the individual particles is becoming harder to see due to the larger densities of particles per given volume. We attempted to provide alternate colors for the plotted particles but Mathematica did not seem to have this capability (AVS was unusable for this large amount of points).

For the last run of 400k on the O2 SG workstation, the rendering in Mathematica was extremely slow (over 2 hours). The cluster had 37949 aggregated particles for the C version and is shown in figure 8.

Once all the DLA clusters were obtained, the fractal analysis for these clusters could be undertaken. The code was modified to output the length of the box and calculate the

Figure 5 - 200K iterations for dla3d.f90

densities per box knowing itís volume and the number of particles enclosed in the box. The length of a given box together with the densities were then plotted in Mathematica to

obtain a log ñ log plot and to obtain the slope of the linear fit of the plotted points. This allowed us to analyse it and to determine if our generated clusters had a fractal structure. From the slope of the density ( versus L, if the cluster is a fractal, then the graph should be a straight line of slope df ñ 2 and ( is proportional to L df ñ 2.

Landau obtains a slope of ñ0.36, which corresponds to a fractal dimension of 1.66. Because random numbers are involved, the graphs we generated will be different, but the

Figure 6 - 300K iterations for dla3d.c

Figure 7 - 300K iterations for dla3d.f90

fractal dimension should be close. Landau notes that the structures for these DLA clusters is multifractal, and to analyse these clusters in detail, a more sophisticated theory and simulation should be used. Due to time constraints, a more sophisticated approach to our analysis wasnít possible.

A number of graphs were obtained and we first analysed the 2D cluster for the 100K iteration shown in Appendix H. Using Landauís box counting approach, the graph shown in figure 10 was obtained. The slope obtained for the graph was ñ 0.52 which corresponds to a fractal dimesion of 1.48. We see that the ìSierpinski gasketî[1] has a dimension between that of a 1D line and a 2D triangle, hence is has a fractional dimension and is a fractal. A number of 2D clusters were generated and the clusters obtained were typical of the one shown in Appendix H. For 50K iterations, a fractional dimension of 1.57 was obtained. For 80K iterations, df was found to be 1.51. From theory df is given to be 1.58496 for the cluster to have a fractal structure as shown in euqation 5.

Figure 8 - 400K iterations for dla3d.f90

Figure 9 - 100K iteration log-log plot for densities versus box length

Having confirmed the code produced clusters of fractal nature, the code was modified to obtain densities using boxes instead of squares for the 2D version. The same procedure is to find df was used but the fractal dimension is calculated using the 3D version of the formula. The slopes were obtained from the log ñlog graph and the fractional dimensions obtained for our 3D clusters are shown in table 1 for DLA3D.C. Similarly the fractional dimensions for the clusters generated using the fortran90 version of our code are shown in table 2.

Iterations df
100K 1.23
200K 1.37
300K 1.51
400K 1.73

Table 1 - Fractional dimensions for given iterations for the 3D cluster runs in C

Iterations df
100K 1.93
200K 1.61
300K 1.74
400K 1.82

Table 2 - Fractional dimensions for given iterations for the 3D cluster runs in Fortran90.

Due to time constraints, the fractal analysis for the multimillion iterations runs on Napier have been left out. It is however safe to assume from the previous results that 5 million and higher iteration runs have produced clusters with fractal structure.

During the Napier runs, it was remarked that using a higher number of processors resulted in the computations taking more time to complete then using less processors and required a substantial larger amount of memory when using more processors.

Discussion

The applications of our 3D DLA cluster growths can have numerous applications in the real world. One example is in explaining planetary formation or the formation of extrasolar planets in other star systems.

Planets in our system are known to have formed from raw materials of the primordial gas nebula from which our Sun was formed. Current accepted theories indicate that stars are born from the gravitational contraction of interstellar coulds of gas and dust. A circumstellar disk around the star results when this gravitational collapse occurs and smaller bodies such as asteroids form from the collisions of dust grains in the circumstellar disk to eventually form planetoids and planets.

Observations of the chemical composition of dense clouds near young stars with circumstellar disks indicate the presence of various elements familiar on Earth. The gaseous component of clouds were found to inlcude hydrogen, helium, nitrogen, oxygen, carbon, sulfur and other compounds while observations of the dust component suggest the presence of silicates, magnetite, carbon compounds such as graphite and water ice all of which are present in our solar system and other planetary systems. The formation of planets, planetesimals and smaller planetary system objects such as small rocks, asteroids etc initially relied on the sticking or coagulation and collision of micrometer sized dust grains thoughout the plane of the circumstellar disk before gravitational accretion by larger objects took place. The growth of these objects from dust grains rely on sticking mechanisms which is largely due to van der Waals bonding and electrostatic forces between the micrometer size dust grains. These forces are relatively independant to the

Figure 10 - Computer simulation of dust grains aggregating into a larger fractal shape object resulting from collisions. This fractal structure results when dust grains in the circumstellar disk stick together to form larger objects. Note the similarity to some of our 3D DLA fractal objects.

composition of the particles. Van der Waals forces alone are enough to produce ~cm size grains in the protoplanetary disk. The mechanisms for grain growth from ~cm to ~m scale objects depends largely on the stickiness, mechanical strength and the fractal dimension of the object, figure 10.

Relative velocities and size differences between the two colliding objects will determine whether they rebound, erode or even fragment. Recent studies have indicated that coagulation for ~cm sized objects occured with sticking probabilities of a few ten percent for particle velocities of a few meters per second. This relatively slow process occurs over a period of 107 years from cloud collapse. The current models indicate that meter sized objects can readily form from similar mechnisms which means there is no obstacle in the formation of ~km size planetesimals. Gravitiational forces between these planetesimals will allow these objects to coagulate to form larger objects such as moons and planets.

The application of our project for planetary formation is obvious. With more work, the physics of the coagulation process of the dust grains in the circumstellar disk can be implemented in our code taking into account the chemical composition of these dust grains and other factors such as mean velocities, temperature etc.

The fractal cluster shown in figure 10 and on the cover of this report are a result of this simulation with the physics implemented which was produced by Levy and Lunine from the University of Arizona. Once they ran the simulation of the cluster growth, they applied a fractal analysis and determined the objectís fractal dimension using our same method and concluded it had a fractal structure. The finding that planets form initially from clusters of dust grains with fractal stuctures is very interesting indeed. To date it is widely accepted that circumstellar disk formation is an integral part of the stellar formation process and there are billions of stars in our galaxy alone and billions of galaxies in our Universe.

Conclusion

This project has investigated diffusion-limited aggregations of fractal structures for 2D and 3D clusters and showed 3D versions of the clusters for various iterations. To confirm their fractal