project004.py 2.42 KB
Newer Older
Litwin Bartosz's avatar
Litwin Bartosz committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
#ising with numba
import numba   
import time
import rich
import argparse
import numpy as np
import random as rn
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from tqdm import tqdm


rich.get_console().clear()

parser = argparse.ArgumentParser(description='ising model simulation.')
parser.add_argument('size', type=int, help='Size of the grid')
parser.add_argument('J', type=float, help='Interaction energy')
parser.add_argument('beta', type=float, help='Inverse temperature')
parser.add_argument('B', type=float, help='External magnetic field')
parser.add_argument('steps', type=int, help='Number of simulation steps')
parser.add_argument('--density', type=float, default=0.5, help='Initial density of spins')
args = parser.parse_args()

def measure_and_statiscs(func):
    def wrapper(*args, **kwargs):
        starttime = time.perf_counter()
        result = func(*args, **kwargs)
        endtime = time.perf_counter()
        print(f"Time measured: {endtime - starttime} seconds")
        return result
    return wrapper


size = args.size
J = args.J
beta = args.beta
B = args.B
steps = args.steps
density = args.density
grid = np.ones((size, size))

@numba.njit
def dE(i,j):
    return 2*J * (grid[i,j] * grid[(i+1)%len(grid),j]+grid[i,j] * grid[(i-1)%len(grid),j]+grid[i,j] * grid[i,(j+1)%len(grid)]+grid[i,j] * grid[i,(j-1)%len(grid)]) + 2*B*grid[i,j]



@numba.njit
def metropolis(grid):
    for i in range(size**2):
        a = rn.randint(0, size - 1)
        b = rn.randint(0, size - 1)
        dE = 2*J * (grid[a,b] * grid[(a+1)%len(grid),b]+grid[a,b] * grid[(a-1)%len(grid),b]+grid[a,b] * grid[a,(b+1)%len(grid)]+grid[a,b] * grid[a,(b-1)%len(grid)]) + 2*B*grid[a,b]

        #dE(a,b)
        if dE < 0 or np.exp(-beta * dE) > rn.random():
            grid[a, b] *= -1


def M():
    return np.sum(grid)/(size*size)


def run():
    starttime = time.perf_counter()
    grid = np.random.choice([1, -1], size=(size, size), p=[density, 1-density])
    fig = plt.figure()
    ims = []
    file = open('magnetyzacja.txt', 'w')
    for i in range(steps):
        metropolis(grid)
        file.write("krok: " + str(i) + " magnetyzacja: " + str(M()) + "\n")
        im = plt.imshow(grid, animated=True)
        ims.append([im])
    file.close()
    ani = animation.ArtistAnimation(fig, ims, interval=50, blit=True, repeat_delay=1000)
    endtime = time.perf_counter()
    print(f"Time measured: {endtime - starttime} seconds")
    plt.show()


run()