Numerical simulation in Python-Invasive percolation
Introduction
Since Python can really do anything, I think there are many ways to use it, but I would like to introduce how I use it. A commonly used package is numpy. It is a companion to numerical calculation. I'm really in trouble without this. The rest is matplotlib and scipy. Then I'm using Tkinter to display the dialog. I think that there are quite a lot of people using this area, so I think it's easy to collect information.
Finally, I would like to create an invasive percolation cluster and put a whole program to calculate its fractal dimension. By default, numpy, scipy, and matplotlib are not included, so you need to include them first.
When I was in a hurry to get into my laptop in the lab
windows7 --Python, SciPy, Install matplotlib (Windows) --Qiita:
It was insanely convenient. Thank you for your help.
When you put it in Ubuntu, you should use apt obediently.
For the model of the program below, please forgive me by copying and pasting my report ...
About invasive percolation
A dynamic process known as invasive percolation is used to model the shape of the oil-water interface that occurs when water is pushed into a porous medium containing oil. The idea is to use water to recover as much oil as possible. In this process, water clusters grow into the oil through the least resistant pathway. Consider an L × 2L grid, and assume that the left edge is initially occupied by water (invaders). The strength of resistance to the invaders is given to each point in the grid with a uniform random number between 0 and 1 and remains fixed during the process of invasion. The closest points to the grid points of all invaders are the surrounding points. At each division time, the surrounding points with the smallest random numbers are occupied by the invaders, and the oil (defender) at that point is replaced. The invader's cluster grows until a pathway connecting the left and right ends of the grid is formed. Periodic boundary conditions are used for the top and bottom edges to minimize the effects of boundaries, and all quantities are measured only in the L × L region of the center of the grid. The main quantities of interest are the percentage of grid points occupied by the invaders and the probability P (r) that grid points with random numbers between r + dr are occupied.
So, now that I've explained the model image, I'd like to touch on the content of the code.
Code description
First, from the basics
The grid points are represented like a matrix using numpy's array (where it says self.lattice).
numpy.random.random([x, y])
Creates an array of random numbers with elements 0 to 1 in a matrix of x rows and y columns. In this model, it is assumed that the value of the random number given to each represents the resistance. Next, as we are doing in the while clause, as described in the comment, the grid points with the minimum resistance value are occupied in order from the points around the occupied grid points.
Where I devised
#Set the coordinates of the grid with the smallest value at the surrounding points to mm
mm = min([(v,k) for k,v in nnsite.items()])[1]
The key and value of the dictionary nnsite are taken out by nnsite.items (), and the key with the minimum value is obtained by comparing the magnitude of the first value of the passed tuple, which is the specification of min.
class TopWindow:
def quit(self):
self.root.destroy()
sys.exit()
def show_window(self, title="title", *args):
self.root = Tk()
self.root.title(title)
frames = []
for i, arg in enumerate(args):
frames.append(Frame(self.root, padx=5, pady=5))
for k, v in arg:
Button(frames[i],text=k,command=v).pack(expand=YES, fill='x')
frames[i].pack(fill='x')
f = Frame(self.root, padx=5, pady=5)
Button(f,text='quit',command=self.quit).pack(expand=YES, fill='x')
f.pack(fill='x')
self.root.mainloop()
Specify the function associated with the button label with the tuple set given to the argument. When calling
run = (('run', pushed), ('save canvas to sample.eps', pr))
run2 = (('calculate D', b4_pushed),)
top.show_window("Invasion Percolation", run, run2)
To do.
I can't explain all the other details, so I will introduce what kind of screen is obtained as a result of the simulation and what kind of execution result is obtained.
In this way, a dialog will be displayed, clicking run will generate a percolation cluster, and calculate_D will give you the graph above. The meaning is omitted.
invasion_percolation.py
#! /usr/bin/env python
# -*- coding:utf-8 -*-
#
# written by ssh0, August 2014.
from Tkinter import *
import sys
import numpy as np
import scipy.optimize as optimize
import matplotlib.pyplot as plt
class Percolation:
def __init__(self, Lx=40, Ly=20):
self.sub = None
self.Lx = Lx
self.Ly = Ly
def perc_cluster(self):
self.lattice = np.random.random([self.Lx+2, self.Ly])
Lx = self.Lx
Ly = self.Ly
#The left end is all occupied sites
self.lattice[:1, :] = 1.5
self.lattice[Lx+1:, :] = 0
if self.sub is None or not self.sub.winfo_exists():
lattice = self.lattice
ne = [(0, -1), (0, 1), (-1, 0), (1, 0)]
#Keep the coordinates of surrounding points in dictionary format
nnsite = {(1, y):lattice[1, y] for y in range(Ly)}
percolate = False
while len(nnsite) != 0 and percolate == False:
#Set the coordinates of the grid with the smallest value at the surrounding points to mm
mm = min([(v,k) for k,v in nnsite.items()])[1]
lattice[mm] = 1
del nnsite[mm]
#List the coordinates of the points around mm in nn(Apply periodic boundary conditions in the y direction)
nn = [(mm[0] + nx, (mm[1] + ny)%Ly) for nx, ny in ne
if lattice[mm[0] + nx, (mm[1] + ny)%Ly] < 1]
#Excludes nn that are already included in nnsite--> nnlist
nnlist = list(set(nn) - set(nnsite.keys()))
#Add new peripheral points to nnsite
for n in nnlist:
nnsite[n] = lattice[n]
#When the occupied grid point is at the right end, it is considered to be percolated.
if mm[0] == Lx:
percolate = True
self.lattice = lattice[1:-1, :]
return self.lattice
def draw_canvas(self, rect, Lx, Ly):
default_size = 640 # default size of canvas
r = int(default_size/(2*Lx))
fig_size_x = 2*r*Lx
fig_size_y = 2*r*Ly
margin = 10
sub = Toplevel()
sub.title('invasion percolation')
self.canvas = Canvas(sub, width=fig_size_x+2*margin,
height=fig_size_y+2*margin)
self.canvas.create_rectangle(margin, margin,
fig_size_x+margin, fig_size_y+margin,
outline='black', fill='white')
self.canvas.pack()
c = self.canvas.create_rectangle
site = np.where(rect==1)
for m, n in zip(site[0], site[1]):
c(2*m*r+margin, 2*n*r+margin,
2*(m+1)*r+margin, 2*(n+1)*r+margin,
outline='black', fill='black')
def get_fractal_dim(self, trial=20, Lmin=20, Lmax=40, Lsample=10):
#An integer value between Lmin and Lmax so that it is evenly spaced when logged as much as possible.
L = np.array([int(i) for i
in np.logspace(np.log10(Lmin), np.log10(Lmax), Lsample)])
M_L = []
for l in L:
self.Lx = l*2
self.Ly = l
m_L = 0
for i in range(trial):
lattice = self.perc_cluster()
#Total number of occupied sites in the central L × L grid
m_L += np.sum(lattice[int(l/2)+1:l+int(l/2),:]==1)
M_L.append(m_L/float(trial))
print "L = %d, M_L = %f" % (l, M_L[-1])
M_L = np.array(M_L)
def fit_func(parameter0, L, M_L):
log = np.log
c1 = parameter0[0]
c2 = parameter0[1]
residual = log(M_L) - c1 - c2*log(L)
return residual
parameter0 = [0.1, 2.0]
result = optimize.leastsq(fit_func, parameter0, args=(L, M_L))
c1 = result[0][0]
D = result[0][1]
print "D = %f" % D
def fitted(L, c1, D):
return np.exp(c1)*(L**D)
fig = plt.figure("Fractal Dimension")
ax = fig.add_subplot(111)
ax.plot(L, M_L, '-o', label=r"$M(L)$")
ax.plot(L, fitted(L, c1, D), label="fit func: D = %f" % D)
ax.set_xlabel(r'$\ln L$', fontsize=16)
ax.set_ylabel(r'$\ln M(L)$', fontsize=16)
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_ymargin(0.05)
fig.tight_layout()
plt.legend(loc='best')
plt.show()
class TopWindow:
def quit(self):
self.root.destroy()
sys.exit()
def show_window(self, title="title", *args):
self.root = Tk()
self.root.title(title)
frames = []
for i, arg in enumerate(args):
frames.append(Frame(self.root, padx=5, pady=5))
for k, v in arg:
Button(frames[i],text=k,command=v).pack(expand=YES, fill='x')
frames[i].pack(fill='x')
f = Frame(self.root, padx=5, pady=5)
Button(f,text='quit',command=self.quit).pack(expand=YES, fill='x')
f.pack(fill='x')
self.root.mainloop()
if __name__ == '__main__':
Lx = 40
Ly = 20
top = TopWindow()
per = Percolation(Lx, Ly)
count = 1
def pr():
global count
d = per.canvas.postscript(file="figure_%d.eps" % count)
print "saved the figure to a eps file"
count += 1
def pushed():
per.perc_cluster()
per.draw_canvas(per.lattice, Lx, Ly)
def b4_pushed():
trial = 100; Lmin = 20; Lmax = 100; Lsample = 10
per.get_fractal_dim(trial, Lmin, Lmax, Lsample)
run = (('run', pushed), ('save canvas to sample.eps', pr))
run2 = (('calculate D', b4_pushed),)
top.show_window("Invasion Percolation", run, run2)