Newer
Older
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
import numpy as np
import sys
from os.path import isfile
from parallel import Parallel
from heat_io import read_restart, read_field
def initialize():
timesteps = 500 # Number of time-steps to evolve system
nx = 1000
ny = 1000
if isfile('HEAT_RESTART.dat'):
field, field0, parallel, iter0 = read_restart()
else:
if len(sys.argv) == 2:
filename = sys.argv[1]
field, field0, parallel = read_field(filename)
elif len(sys.argv) == 3:
filename = sys.argv[1]
timesteps = int(sys.argv[2])
field, field0, parallel = read_field(filename)
elif len(sys.argv) == 4:
nx = int(sys.argv[1])
ny = int(sys.argv[2])
timesteps = int(sys.argv[3])
field, field0, parallel = generate_field(nx, ny)
else:
field, field0, parallel = generate_field(nx, ny)
iter0 = 0
return field, field0, parallel, iter0, timesteps
def generate_field(nx, ny):
parallel = Parallel(nx, ny)
nx_local = nx / parallel.dims[0]
ny_local = ny / parallel.dims[1]
field = np.zeros((nx_local + 2, ny_local + 2), dtype=float)
# Square of the disk radius
radius2 = (nx / 6.0)**2
X, Y = np.meshgrid(np.arange(nx_local + 2), np.arange(ny_local + 2),
indexing='ij')
ds2 = (X + parallel.coords[0] * nx_local - nx / 2.0 )**2 + \
(Y + parallel.coords[1] * ny_local - ny / 2.0 )**2
field[:] = np.where(ds2 < radius2, 5.0, 65.0)
# Boundaries
if parallel.coords[0] == 0:
field[0,:] = 85.0
if parallel.coords[0] == parallel.dims[0] - 1:
field[-1,:] = 5.0
if parallel.coords[1] == 0:
field[:,0] = 20.0
if parallel.coords[1] == parallel.dims[1] - 1:
field[:,-1] = 70.0
field0 = field.copy()
return field, field0, parallel