Skip to content

Commit cb1970d

Browse files
author
Kyle Schau
committed
Working.
1 parent 274e7d7 commit cb1970d

File tree

3 files changed

+228
-208
lines changed

3 files changed

+228
-208
lines changed

utilities/animateAlphas.py

Lines changed: 101 additions & 91 deletions
Original file line numberDiff line numberDiff line change
@@ -25,12 +25,13 @@
2525
import numpy as np
2626
import peregrinepy as pg
2727
import yaml
28-
from mpl.animation import FuncAnimation
28+
from matplotlib.animation import FuncAnimation
2929

30-
inputFile = sys.argv[1]
3130

31+
inputFile = sys.argv[1]
32+
bcFam = inputFile.split(".")[0]
3233
# Use the input file name as the alpha directory
33-
outputDir = inputFile.split(".")[0] + "_alphas"
34+
outputDir = bcFam + "_alphas"
3435
if not os.path.exists(outputDir):
3536
raise IOError(f"Cannot find your alpha files in the usual place, {outputDir}.")
3637

@@ -39,49 +40,51 @@
3940
###############################################################################
4041
with open(inputFile, "r") as f:
4142
seminp = yaml.load(f, Loader=yaml.FullLoader)
42-
43-
###############################################################################
44-
# Read in patches.txt
45-
###############################################################################
46-
# Only the zeroth rank reads in patches.txt
47-
patchNum = []
48-
blockNum = []
49-
with open("patches.txt", "r") as f:
50-
for line in [i for i in f.readlines() if not i.startswith("#")]:
51-
li = line.strip().split(",")
52-
patchNum.append(int(li[0]))
53-
blockNum.append(int(li[1]))
54-
55-
npatches = len(patchNum)
56-
43+
nframes = seminp["nframes"]
5744
###############################################################################
5845
# PEREGRINE stuff (read in grid and make the patches)
5946
###############################################################################
6047
# Only rank ones reads in the grid
61-
nblks = len([f for f in os.listdir(seminp["gridPath"]) if f.startswith("g.")])
48+
nblks = len(
49+
[
50+
f
51+
for f in os.listdir(seminp["gridPath"])
52+
if f.startswith("gv.") and f.endswith(".h5")
53+
]
54+
)
6255
if nblks == 0:
6356
nblks = len([f for f in os.listdir(seminp["gridPath"]) if f == "grid"])
6457
if nblks == 0:
6558
raise ValueError(f'Cant find any grid files in {seminp["gridPath"]}')
6659

6760
mb = pg.multiBlock.grid(nblks)
6861
pg.readers.readGrid(mb, seminp["gridPath"])
62+
pg.readers.readConnectivity(mb, seminp["connPath"])
63+
inletBlocks = []
64+
inletFaces = []
65+
for blk in mb:
66+
for face in blk.faces:
67+
if face.bcFam == bcFam:
68+
inletBlocks.append(blk.nblki)
69+
inletFaces.append(face.nface)
70+
print(f"\nFound {len(inletBlocks)} blocks for this inlet.")
71+
6972
# Flip the grid if it is oriented upside down in the true grid
7073
if seminp["flipdom"]:
71-
for bn in blockNum:
72-
blk = mb[bn - 1]
74+
for bn in inletBlocks:
75+
blk = mb.getBlock(bn)
7376
blk.array["y"] = -blk.array["y"]
7477
blk.array["z"] = -blk.array["z"]
7578
# Determinw the extents the grid needs to be shifted
7679
ymin = np.inf
7780
zmin = np.inf
78-
for bn in blockNum:
79-
blk = mb[bn - 1]
81+
for bn in inletBlocks:
82+
blk = mb.getBlock(bn)
8083
ymin = min(ymin, blk.array["y"][0, :, :].min())
8184
zmin = min(zmin, blk.array["z"][0, :, :].min())
8285
# Now shift the grid to match the domain (0,0) starting point
83-
for bn in blockNum:
84-
blk = mb[bn - 1]
86+
for bn in inletBlocks:
87+
blk = mb.getBlock(bn)
8588
blk.array["y"] = blk.array["y"] - ymin
8689
blk.array["z"] = blk.array["z"] - zmin
8790
# Compute the locations of face centers
@@ -94,47 +97,26 @@
9497
y = np.array([])
9598
z = np.array([])
9699

97-
ny = []
98-
nz = []
99-
for bn in blockNum:
100-
blk = mb[bn - 1]
101-
y = np.concatenate((y, blk.array["iyu"][0, :, :].ravel()))
102-
z = np.concatenate((z, blk.array["izu"][0, :, :].ravel()))
103-
ny.append(blk.nj)
104-
nz.append(blk.nk)
100+
for bn, fn in zip(inletBlocks, inletFaces):
101+
blk = mb.getBlock(bn)
102+
face = blk.getFace(fn)
103+
if fn == 1:
104+
s1_ = np.s_[0, :, :]
105+
elif fn == 2:
106+
s1_ = np.s_[-1, :, :]
107+
elif fn == 3:
108+
s1_ = np.s_[:, 0, :]
109+
elif fn == 4:
110+
s1_ = np.s_[:, -1, :]
111+
elif fn == 5:
112+
s1_ = np.s_[:, :, 0]
113+
elif fn == 6:
114+
s1_ = np.s_[:, :, -1]
115+
y = np.concatenate((y, blk.array["iyc"][s1_].ravel()))
116+
z = np.concatenate((z, blk.array["izc"][s1_].ravel()))
105117

106118
tri = mpl.tri.Triangulation(z, y)
107119

108-
###############################################################################
109-
# Function to read in frames
110-
###############################################################################
111-
112-
113-
def getFrame(frame, comp):
114-
alphas = np.array([])
115-
for i, pn in enumerate(patchNum):
116-
fileName = outputDir + "/alphas_{:06d}_{:04d}".format(pn, frame)
117-
with open(fileName, "rb") as f:
118-
# The coefficients are in "reverse" order, i.e. the last
119-
# of the four coefficients is the constant value at the
120-
# t-t[i] frame location. So we just plot that.
121-
if comp == "u":
122-
a = np.load(f)[-1, :, :]
123-
elif comp == "v":
124-
_ = np.load(f)
125-
a = np.load(f)[-1, :, :]
126-
elif comp == "w":
127-
_ = np.load(f)
128-
_ = np.load(f)
129-
a = np.load(f)[-1, :, :]
130-
else:
131-
raise ValueError(f"Unknown component {comp}")
132-
133-
alphas = np.concatenate((alphas, a.ravel()))
134-
135-
return alphas
136-
137-
138120
###############################################################################
139121
# Make plots
140122
###############################################################################
@@ -146,32 +128,50 @@ def getFrame(frame, comp):
146128
plt.rcParams["font.size"] = 12
147129
plt.rcParams["image.cmap"] = "autumn"
148130

149-
# animation function
150-
151-
152-
def animate(i, Tri, comp):
131+
alphasU = np.array([]).reshape(nframes - 1, 0)
132+
alphasV = np.array([]).reshape(nframes - 1, 0)
133+
alphasW = np.array([]).reshape(nframes - 1, 0)
134+
for bn, fn in zip(inletBlocks, inletFaces):
135+
fileName = outputDir + "/alphas_{}_{}.npy".format(bn, fn)
136+
temp = np.array([])
137+
with open(fileName, "rb") as f:
138+
# The coefficients are in "reverse" order, i.e. the last
139+
# of the four coefficients is the constant value at the
140+
# t-t[i] frame location. So we just plot that.
141+
au = np.load(f)[-1, :, :, :]
142+
av = np.load(f)[-1, :, :, :]
143+
aw = np.load(f)[-1, :, :, :]
144+
shape = au.shape
145+
au = au.ravel().reshape((nframes - 1, np.prod(au.shape[1::])))
146+
av = av.ravel().reshape((nframes - 1, np.prod(av.shape[1::])))
147+
aw = aw.ravel().reshape((nframes - 1, np.prod(aw.shape[1::])))
148+
alphasU = np.hstack((alphasU, au))
149+
alphasV = np.hstack((alphasV, av))
150+
alphasW = np.hstack((alphasW, aw))
151+
152+
153+
def animate(i, Tri, alphas, comp):
153154
global tcf
154155
for c in tcf.collections:
155156
c.remove() # removes only the contours, leaves the rest intact
156-
alphas = getFrame(i + 1, comp)
157-
tcf = ax.tricontourf(Tri, alphas, levels=levels)
157+
tcf = ax.tricontourf(Tri, alphas[i], levels=levels)
158158
ax.set_title(f"Frame {i+1}")
159159
return tcf
160160

161161

162162
###############################################################################
163163
# Film u
164164
###############################################################################
165+
plt.rcParams["image.cmap"] = "plasma"
165166
print("Animating U")
166167
fig, ax = plt.subplots()
167-
ax.set_xlabel(r"$z/L_{ref}$")
168-
ax.set_ylabel(r"$y/L_{ref}$")
168+
ax.set_xlabel(r"$z$")
169+
ax.set_ylabel(r"$y$")
169170
ax.set_aspect("equal")
170171

171-
alphas = getFrame(1, "u")
172-
levels = np.linspace(0, alphas.max(), 100)
173-
ticks = np.linspace(0, alphas.max(), 11)
174-
tcf = ax.tricontourf(tri, np.zeros(alphas.shape), levels=levels)
172+
levels = np.linspace(0, alphasU.max(), 100)
173+
ticks = np.linspace(0, alphasU.max(), 11)
174+
tcf = ax.tricontourf(tri, np.zeros(alphasU[0].shape), levels=levels)
175175
ratio = y.max() / z.max()
176176
cb = plt.colorbar(
177177
tcf,
@@ -183,7 +183,11 @@ def animate(i, Tri, comp):
183183

184184

185185
anim = FuncAnimation(
186-
fig, animate, frames=seminp["nframes"] - 1, fargs=(tri, "u"), repeat=False
186+
fig,
187+
animate,
188+
frames=seminp["nframes"] - 1,
189+
fargs=(tri, alphasU, "u"),
190+
repeat=False,
187191
)
188192
try:
189193
anim.save("U.mp4", writer=mpl.animation.FFMpegWriter(fps=10))
@@ -199,26 +203,27 @@ def animate(i, Tri, comp):
199203
print("Animating v")
200204

201205
fig, ax = plt.subplots()
202-
ax.set_xlabel(r"$z/L_{ref}$")
203-
ax.set_ylabel(r"$y/L_{ref}$")
206+
ax.set_xlabel(r"$z$")
207+
ax.set_ylabel(r"$y$")
204208
ax.set_aspect("equal")
205209

206-
alphas = getFrame(1, "v")
207-
symm = max(alphas.max(), abs(alphas.min()))
210+
symm = max(alphasV.max(), abs(alphasV.min()))
208211
levels = np.linspace(-symm, symm, 100)
209212
ticks = np.linspace(-symm, symm, 11)
213+
tcf = ax.tricontourf(tri, np.zeros(alphasV[0].shape), levels=levels)
210214
cb = plt.colorbar(
211-
tcf, ticks=ticks, fraction=0.046 * ratio, pad=0.05, label=r"$v^{\prime}$"
215+
tcf,
216+
ticks=ticks,
217+
fraction=0.046 * ratio,
218+
pad=0.05,
219+
label=r"$v^{\prime}$",
212220
)
213221

214-
anim = mpl.animation.FuncAnimation(
222+
anim = FuncAnimation(
215223
fig,
216224
animate,
217225
frames=seminp["nframes"] - 1,
218-
fargs=(
219-
tri,
220-
"v",
221-
),
226+
fargs=(tri, alphasV, "v"),
222227
repeat=False,
223228
)
224229
try:
@@ -234,17 +239,21 @@ def animate(i, Tri, comp):
234239
print("Animating w")
235240

236241
fig, ax = plt.subplots()
237-
ax.set_xlabel(r"$z/L_{ref}$")
238-
ax.set_ylabel(r"$y/L_{ref}$")
242+
ax.set_xlabel(r"$z$")
243+
ax.set_ylabel(r"$y$")
239244
ax.set_aspect("equal")
240245

241-
alphas = getFrame(1, "w")
242-
symm = max(alphas.max(), abs(alphas.min()))
246+
symm = max(alphasW.max(), abs(alphasW.min()))
243247
symm = float(f"{symm:.2e}")
244248
levels = np.linspace(-symm, symm, 100)
245249
ticks = np.linspace(-symm, symm, 11)
250+
tcf = ax.tricontourf(tri, np.zeros(alphasW[0].shape), levels=levels)
246251
cb = plt.colorbar(
247-
tcf, ticks=ticks, fraction=0.046 * ratio, pad=0.05, label=r"$w^{\prime}$"
252+
tcf,
253+
ticks=ticks,
254+
fraction=0.046 * ratio,
255+
pad=0.05,
256+
label=r"$w^{\prime}$",
248257
)
249258

250259
anim = mpl.animation.FuncAnimation(
@@ -253,6 +262,7 @@ def animate(i, Tri, comp):
253262
frames=seminp["nframes"] - 1,
254263
fargs=(
255264
tri,
265+
alphasW,
256266
"w",
257267
),
258268
repeat=False,

0 commit comments

Comments
 (0)