Homework: DLA’ing with fractals

• Topic: random walks, dynamics, growth, fractals

• Template: A2_dla.py

• Further reading:

  • https://en.wikipedia.org/wiki/Diffusion-limited_aggregation

  • https://en.wikipedia.org/wiki/Fractal_dimension

  • http://paulbourke.net/fractals/dla/

Theory

Diffusion limied aggregation (DLA) is a growth model operating on the following principle: particles in, say, a liquid move randomly. These particles can accumulate together to form a cluster, if they happen to meet each other. However, as the cluster grows, it is most likely that the diffusing particles attach to the outer edge of the cluster, and so the cluster gets a particular fractal-like shape.

In this assignment we simulate the formation of a discrete 2D DLA cluster. That is, we simulate randomly diffusing particles on a lattice plane with a cluster in the center. Once a diffusing particle meets the cluster, it attaches to it and cannot move any longer. Then, another particle is created and allowed to diffuse, and so on.

DLA clusters have the same relative roughness at different length scales, meaning that a cluster of \(10^9\) particles looks pretty much the same as a cluster of \(10^3\) particles, only bigger. This is a property of fractals, and it is quantified by the fractal dimension.

The basic idea of a fractal dimension is this: If one draws a line of width \(d\) (a 1D object) from a given point and then draws circle of radius \(R\) around the starting point, the area inside the circle covered by the line is \(A_\text{line} = d R\). On the other hand if one draws a triangle (a 2D object) with one vertex at the given point, this triangle covers the area \(A_\text{triangle} = \frac{\phi}{2}R^2\), where \(\phi\) is the angle of the triangle vertex. The total area of the circle is of course \(A_\text{circle} = \pi R^2\).

Now, consider the logarithm of the covered area. We get for the line

\[\ln A_\text{line} = \ln d R = \ln d + \ln R\]

and for the triangle

\[\ln A_\text{triangle} = \ln \frac{\phi}{2} R^2 = \ln \frac{\phi}{2} + 2 \ln R.\]

In general, for a traditional geometric object, the logarithm of the area (or volume) grows linearly as a function of the logarithm of the radius with a slope equal to the dimensionality of the object, \(D\),

\[\ln A = C + D \ln R.\]

For a line, \(D = 1\), and for a planar shape such as a triangle, \(D = 2\). The above relationship also holds for a fractal, but the dimension of a fractal can be non-integer.

Furthermore, this result lets us determine the dimensionality of a fractal. If we plot \(\ln A\), or in the case of our simulation, \(\ln n\) where \(n(R)\) is the number of particles inside a circle of radius \(R\), as a function of \(\ln R\), we get a straight line whose slope is the fractal dimension \(D\).

Task

To complete this assignment, you need to complete the following tasks, submit your code as a plain .py file, and answer the questions. In the computing tasks, you must implement the incomplete functions so that they operate exactly as described in the documentation.

• Computing task:

  1. Implement a function to calculate the radius of the DLA cluster.

  2. Implement a function to simulate the diffusion of a single particle.

  3. Implement a function to simulate the diffusion of several particles and thus the formation of the DLA cluster.

  4. Implement a function to calculate the fractal dimension of the DLA cluster. (It will be tested for clusters of 5000-30000 particles.)

• Questions to answer:

  1. How does your program find \(D\)? Give a detailed explanation. Include a picture of your cluster and the plot you used in the analysis.

  2. What is the fractal dimension of a DLA cluster according to your results?

  3. Investigate the linked material. What should the fractal dimension be? (Cite your source!) Is your result close to the literature value?

A2_dla.py

A2_dla.animate(history)

Animate the simulation.

Parameters:

history (list) – list of systems at different times

A2_dla.calculate_cluster_radius(grid, center)

Loops over all the sites in the grid and calculates the distance from the center coordinates. Returns the maximum distance found for occupied sites.

That is, the function calculates the radius of the smallest circle that is centered at the given center coordinates and which overlaps the entire cluster.

Note

This function is incomplete!

Parameters:

• grid (array) – the grid where values of 1 represent the cluster

• center (array) – cluster center coordinates [x0, y0]

Returns:

the calculated cluster radius

Return type:

float

A2_dla.calculate_fractal_dimension(grid, center, printout=True)

Calculate the fractal dimenstion of the cluster and optionally display the plot from which the dimension can be determined.

In \(d\) dimensions, the size \(s\) of an object increases as \(s = cr^d\) as the linear size \(r\) of the object grows, where \(c\) is some constant. For instance, the area of a 2D circle is \(A = \pi r^2\) and the volume of a 3D ball is \(V = \frac{4}{3} \pi r^3\).

Taking the logarithm, we get \(\ln(A) = \ln(\pi) + 2 \ln(r)\) and \(\ln(V) = \ln(\frac{4}{3}\pi) + 3 \ln(r)\) and in general \(\ln(s) = \ln(c) + d \ln(r)\). This means we can find the dimensionality of an object by looking at how its size grows as its linear length grows on a logarithmic scale.

Although ideal mathematical shapes tend to have integer dimensionality, real objects are usually fractals at least on some lenght scales. This means their dimensionality can be non-integer.

The fractal dimension of the cluster is calculated as follows:

• Separate circles of radius R = 1,2,3,… around the cluster center.

• Calculate the number of cluster particles, N, (grid points with value 1) in each circle (using count_particles_in_circle()).

• Optionally: plot the N(R) data on a log-log scale.

• Estimate the region where the data falls on a line.

• Fit a linear function to the data, in this region. The slope of the fitted line is the fractal dimension.

Note

This function is incomplete!

Parameters:

• grid (array) – an array of points where the value of 1 denotes a point that belongs in the cluster

• center (array) – the coordinates of the initial seed of the cluster

• printout (bool) – if True, results are shown

Returns:

the calculated estimate for fractal dimension

Return type:

float

A2_dla.calculate_particles_in_cluster(grid)

Loops over all the sites in the grid and counts the number of sites with value 1.

Parameters:

grid (array) – the grid where values of 1 represent the cluster

Returns:

number of particles

Return type:

int

A2_dla.count_particles_in_circle(center, R_max, grid)

Count the number of points in the grid with the value of 1 for which the distance from center is at most R_max.

That is, count how many of the grid points in a circle with radius R_max belong to the cluster.

Parameters:

• center (array) – cluster center coordinates [x0, y0]

• R_max (float) – maximum distance from center

• grid (array) – the grid to be analyzed

Returns:

the number of gridpoints = 1 in the circle

Return type:

int

A2_dla.count_particles_in_square(center, boxsize, grid)

Count the number of points in the grid with the value of 1 for which the x and y indices differ from the coordinates [x0 ,y0] given in ‘center’ by at most boxsize.

That is, count how many of the grid points in [x0-boxsize, x0+boxsize] x [y0-boxsize, y0+boxsize] are 1.

Parameters:

• center (array) – cluster center coordinates [x0, y0]

• boxsize (int) – maximum distance from center

• grid (array) – the grid to be analyzed

Returns:

the number of gridpoints = 1 in the box

Return type:

int

A2_dla.create_grid(n)

Creates a grid as n x n numpy array.

The grid contains 0’s except the center point is 1. This represents the initial seed for the cluster.

Parameters:

n (int) – size of the grid

Returns:

grid, center coordinates [x0, y0]

Return type:

array, array

A2_dla.diffusion_simulation(grid, center, particles, recording_interval=100)

Run a full diffusion simulation on the given grid.

The grid must be an array of zeros and ones. It represents the space in which the simulation takes place as a discrete lattice model. That is, each particle to be simulated can only have integer coordinates. Each zero in the grid represents empty space. Each one in the grid represents a position filled by the cluster.

The simulation is carried out as follows:

• Particles are released, one-by-one, outside the cluster.

• A particle diffuses until they end up at a position next to the existing cluster via simulate_single_particle().

• When this happens, the particle attaches to the cluster. This is represented by changing the value of the grid to 1 at that position.

• The simulation ends after the given number of particles have joined the cluster.

Note

This function is incomplete!

Parameters:

• grid (array) – a square grid of integer values

• center (array) – center coordinates [x0, y0]

• particles (int) – number of particles to simulate

• recording_interval (int) – record the state of the system after this many particles

Returns:

list of grids showing a time series of the evolution of the system

Return type:

list

A2_dla.draw(frame, history)

Draws the system for animation.

Parameters:

• frame (int) – index of the frame to draw

• history (list) – list of systems at different times

A2_dla.linearfit(x, a, b)

Calculates the linear function \(ax+b\) and returns the result.

Parameters:

• x (float) – the variable

• a (float) – the slope

• b (float) – the constant term

Returns:

the result

Return type:

float

A2_dla.main(filename=None, size=300, n_particles=0)

The main program.

Runs a simulation, shows the resulting cluster and calculates its fractal dimension. The end result is written in a file using write_cluster_file().

Can either start from scratch (a cluster of 1 particle) or read in the result from a previous simulation and use that as the starting point.

If a file is read successfully, the size parameter does nothing. If no file is given or an error occurs while reading, a new size x size simulation grid is created.

Parameters:

• filename (str) – name of the file to read for the starting configuration

• size (int) – the size of the simulation space, if starting from scratch

• n_particles (int) – the number of particles to simulate

A2_dla.mind_the_gap(x, y, center, cluster_radius, margin=10)

Check the distance between the cluster and the particle and decrease it if the distance is too large.

If the particle is far from the cluster, we want to move it a bit closer. This saves simulation time by making sure the particle cannot wander aimlessly away from the cluster. That could lead to very long random walks.

In practice, the function calculates the distance betweeen the given (x,y) coordinates and the center coordinates. If this distance is larger than cluster_radius + margin, the function calculates a new point which is in the same direction from the center as the point (x,y) but 0.5*margin closer.

If the original point was not too far, the function returns the original coordinates. If it was, the function returns the new coordinates rounded to integers values.

Parameters:

• x (int) – the x coordinate

• y (int) – the y coordinate

• center (array) – cluster center coordinates [x0, y0]

• cluster_radius (float) – the maximum radius of the cluster

• margin (float) – the allowed distance in addition to cluster radius

Returns:

new x, new y

Return type:

int, int

A2_dla.neighbor_is_in_cluster(x, y, grid)

Check if any of the neighbors of the given point are already a part of the cluster.

The neighbors of the point (x,y) are the points (x-1, y), (x+1, y), (x, y-1) and (x, y+1).

A point is part of the cluster if that point has the value of 1 in the given array ‘grid’.

Parameters:

• x (int) – the x coordinate

• y (int) – the y coordinate

• grid (array) – the grid to be checked

Returns:

True if any neighbor is in the cluster

Return type:

bool

A2_dla.pick_random_point(center, distance)

Choose a random point at the given distance from the center point.

Distance should be a real number and center should be an array of two real numbers, [x, y]. The result is rounded to integer coordinates.

Parameters:

• center (list) – center point coordinates [x0, y0]

• distance (float) – distance from center

Returns:

integer coordinates as an array [x, y]

Return type:

array

A2_dla.print_progress(step, total)

Prints a progress bar.

Parameters:

• step (int) – progress counter

• total (int) – counter at completion

A2_dla.random_step(x, y)

Take a random step of length 1 from the given coordinates (x,y) and return the new coordinates.

Parameters:

• x (int) – the x coordinate

• y (int) – the y coordinate

Returns:

the new coordinates, [x,y]

Return type:

list

A2_dla.read_grid_from_file(filename)

Read a cluster from a file.

The first line of the file should contain the coordinates of the initial seed of the cluster. This should be followed by the coordinates of all other points in the cluster, line by line:

x0 y0
x1 y1
x2 y2
...

The function returns a grid describing the cluster and the coordinates of the initial seed.

Parameters:

filename (str) – file to be read

Returns:

grid, cluster center coordinates [x0, y0]

Return type:

array, array

A2_dla.show_cluster(grid, cut_corners=False)

Draws the system as a black and white pixel image and saves it to cluster.png.

Parameters:

• grid (array) – grid where the cluster is stored as values of 1

• cut_corners (bool) – If True, do not show the entire grid. Instead only draw the area containing the cluster.

A2_dla.simulate_single_particle(grid, center, cluster_radius, start)

Run a diffusion simulation for a single particle on a grid.

• Place the particle at the given position.

• Let the particle perform a random walk on the grid via random_step().

• If at any time the particle wanders too far away from the cluster, move it closer to save simulation time.

• Let the random walk continue until the walker is next to a grid point with the value 1.

• Change the value of the grid point at the last position of the walker to 1. Physically, grid represents a cluster and this action represents a particle joining the cluster.

• Check if the new particle increases the maximum radius of the cluster. The new radius is returned in the end.

Note

This function is incomplete!

Parameters:

• grid (array) – the grid representing the cluster

• center (list) – the center coordinates of the cluster, [x0, y0]

• cluster_radius (float) – maximum distance between the center and all other points in the cluster

• start (list) – the starting coordinates of the particle, [x, y]

Returns:

the new radius

Return type:

float

A2_dla.write_cluster_file(grid, center)

Write the cluster in a datafile which can be later read in using read_grid_from_file().

Parameters:

• grid (array) – the grid where values of 1 represent the cluster

• center (array) – cluster center coordinates [x0, y0]