-
Notifications
You must be signed in to change notification settings - Fork 0
/
oscillator.py
114 lines (80 loc) · 2.71 KB
/
oscillator.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
"""
Simulate flow through any arbitrary structure: in this case, a fluidic oscillator.
ARGUMENTS
--lat_x, --lat_y : int
overrides default x- and y-dimensions of the entire lattice
--grid_x, --grid_y : int
overrides default x- and y-dimensions of the cartesian processor arrangement. If specified, their product must match the number of processors in use.
INPUTS
oscillator.bmp
monochrome bitmap representation of oscillator
OUTPUTS
oscillator.pkl.gz
Compressed results file. Used for plotting.
@author: AlexR
"""
#%% IMPORTS
import numpy as np
from mpi4py import MPI
import PIL
from PIL import Image
import _pickle as pickle
from lattice import Lattice
from utils import fetch_dim_args, pickle_save
#%% SET PARAMETERS
[lat_x, lat_y], grid_dims = fetch_dim_args(lat_default=[400, 300])
interval = 1
omega = 1.0
inflow = 0.01
outfile = 'oscillator.pkl.gz'
t_recordpoints = [
100, 200, 300, 400, 600, 800, 1000, 2000, 3000
]
#%% SETUP
# use the figure we cut out of the exercise sheet (a monochrome bitmap)
im = Image.open("oscillator.bmp")
# resize it to desired lattice dimensions
walls = np.asarray(im.resize([lat_x, lat_y])).T
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
size = comm.Get_size()
# set up the lattice
lat = Lattice([lat_x, lat_y], grid_dims=grid_dims)
# set walls directly (this won't work in parallel) # lat.walls = walls
#cut up oscillator pattern into per-node pieces
lat.walls = walls[lat.cell_starts[rank, 0]:lat.cell_starts[rank, 0] +
lat.cell_dims[rank, 0], lat.cell_starts[rank, 1]:lat.
cell_starts[rank, 1] + lat.cell_dims[rank, 1]]
# find the gap in the left-hand side
inflow_y = np.logical_not(lat.walls[0, :])
# find the gap in the right-hand side
outflow_y = np.logical_not(lat.walls[-1, :])
if rank == 0:
lat.print_info()
# set up container for data used in the streamplot
flow_hist = np.empty([len(t_recordpoints), lat_x, lat_y, 2])
max_timesteps = max(t_recordpoints)
#%% SIMULATION
for t in range(max_timesteps + 1):
lat.halo_copy()
lat.stream()
# get velocity from our cell
u = lat.u()
# prescribe inflow
if lat.cart.coords[0] == 0:
u[0, inflow_y, 0] = inflow
u[0, inflow_y, 1] = 0
lat.collide(omega=omega, u=u)
# record data
if t in t_recordpoints:
u_snapshot = lat.gather(lat.u())
# necessary for streamplots to work later
u_snapshot[walls] = 0
if rank == 0:
flow_hist[t_recordpoints.index(t)] = u_snapshot
#%% SAVE RESULTS
if rank == 0:
d = dict(
((k, eval(k))
for k in ['lat_x', 'lat_y', 'omega', 'inflow', 't_recordpoints', 'flow_hist', 'walls']))
pickle_save(outfile, d)