-
Notifications
You must be signed in to change notification settings - Fork 0
/
heat simulation.py
122 lines (96 loc) · 3.66 KB
/
heat simulation.py
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
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
# -*- coding: utf-8 -*-
"""
Created on Tue Apr 14 00:09:59 2020
@author: user
"""
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
import copy
import Grad_Div_Curl_Laplace as gdcl
n = 64
map_size = 10
unit_size = map_size/n
GDCL = gdcl.__GDCL__(n, map_size)
time = 100
dt = 0.1
x = np.linspace(-map_size/2, map_size/2, n)
y = np.linspace(-map_size/2, map_size/2, n)
X,Y = np.meshgrid(x, y)
T = [[0 for i in range(n)] for j in range(n)]
for i in range(n):
for j in range(n):
if (i-32)**2+(j-32)**2 < 100:
T[i][j] = 100
if (i-32)**2+(j-63)**2 < 400:
T[i][j] = None
if (i-32)**2+j**2 < 400:
T[i][j] = None
e = [[0 for i in range(n)] for j in range(n)]
for i in range(n):
for j in range(n):
if (i-32)**2+(j-32)**2 < 100:
e[i][j] = 100
Alpha = [[0.05 for i in range(n)] for j in range(n)]
Kapha = [[1 for i in range(n)] for j in range(n)]
all_T_list = []
for t in range(int(time/dt)):
print(t/int(time/dt)*100)
all_T_list.append(copy.deepcopy(T))
L_T = GDCL.Laplacian(T)
U = [[0 for i in range(n)] for j in range(n)]
for i in range(n):
for j in range(n):
try:
T[i][j] += Alpha[i][j]*L_T[i][j]*dt + Kapha[i][j]*e[i][j]*dt
except:
T[i][j] = None
U = []
T_list = []
for t in range(time):
for i in range(n):
for j in range(n):
if all_T_list[t*int(1/dt)][i][j] != None:
U.append(all_T_list[t*int(1/dt)][i][j])
else:
all_T_list[t*int(1/dt)][i][j] = np.nan
T_list.append(all_T_list[t*int(1/dt)])
level = np.linspace(min(U), max(U), 50)
fig = plt.figure(figsize=(10, 8))
ax = plt.axes(xlim=(-map_size/2, map_size/2), ylim=(-map_size/2, map_size/2))
ax.set_title('Heat Simulation 0.0s', fontsize=20)
cp = ax.contourf(X, Y, all_T_list[0], 8, levels = level, cmap = 'hot')
ax.contour(X, Y, all_T_list[0], 8, colors='black',levels = level, linewidth=.5)
fig.colorbar(cp)
def animate(i):
ax.clear()
ax.set_title('Heat Simulation %.1fs'%i, fontsize=20)
ax.contourf(X, Y, all_T_list[i*int(1/dt)], 8, levels = level, cmap = 'hot')
C=ax.contour(X, Y, all_T_list[i*int(1/dt)], 8, colors='black',levels = level, linewidth=.5)
# ax2.clabel(C, inline=1, fontsize=10)
#call the animator. blit=True means only re-draw the parts that have changed.
anim = animation.FuncAnimation(fig, animate,
frames=time, interval=10,
blit=False)
anim.save('dynamic_images3.gif',fps=60)
plt.show()
T_list = np.array(T_list)
from mpl_toolkits.mplot3d import Axes3D
fig2 = plt.figure(figsize=(10, 10))
ax2 = Axes3D(fig2)
ax2.set_title('Heat Simulation 0.0s', fontsize=30)
ax2.plot_surface(X, Y, T_list[0], rstride=1, cstride=1, vmin=min(U), vmax=max(U), cmap='hot')
ax2.contourf(X, Y, T_list[0], zdir='z', offset=min(U), levels = level, cmap='hot')
ax2.set_zlim(-min(U), max(U))
def animate2(i):
ax2.clear()
ax2.set_title('Heat Simulation %.1fs'%i, fontsize=30)
ax2.plot_surface(X, Y, T_list[i], rstride=1, cstride=1, vmin=min(U), vmax=max(U), cmap='hot')
ax2.contourf(X, Y, T_list[i], zdir='z', offset=min(U), levels = level, cmap = 'hot')
ax2.set_zlim(-min(U), max(U))
#call the animator. blit=True means only re-draw the parts that have changed.
anim2 = animation.FuncAnimation(fig2, animate2,
frames=time, interval=10,
blit=False)
anim2.save('dynamic_images4.gif',fps=60)
plt.show()