# Author: Vatsal Sanjay
# vatsalsanjay@gmail.com
# Physics of Fluids
import numpy as np
import os
import subprocess as sp
import matplotlib.pyplot as plt
from matplotlib import rc
import matplotlib
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
from matplotlib.collections import LineCollection
from matplotlib.ticker import StrMethodFormatter
import sys
matplotlib.rcParams['text.usetex'] = True
matplotlib.rcParams['text.latex.preamble'] = [r'']
def gettingFacets(filename, tracer):
print('Getting facets values')
if tracer == 1:
exe = ["./getFacet1", filename]
else:
exe = ["./getFacet2", filename]
p = sp.Popen(exe, stdout=sp.PIPE, stderr=sp.PIPE)
stdout, stderr = p.communicate()
temp1 = stderr.decode("utf-8")
temp2 = temp1.split("\n")
segs = []
skip = False
if (len(temp2) > 1e1):
for n1 in range(len(temp2)):
temp3 = temp2[n1].split(" ")
if temp3 == ['']:
skip = False
pass
else:
if not skip:
temp4 = temp2[n1+1].split(" ")
r1, z1 = np.array([float(temp3[1]), float(temp3[0])])
r2, z2 = np.array([float(temp4[1]), float(temp4[0])])
segs.append(((r1, z1),(r2,z2)))
segs.append(((-r1, z1),(-r2,z2)))
skip = True
print('Got facets values for tracer')
return segs
def gettingfield(filename):
print('Field values')
exe = ["./getData_InsideDrop", filename, str(zmin), str(0), str(zmax), str(rmax), str(nz), str(Ohd), str(Ohf), str(hf)]
p = sp.Popen(exe, stdout=sp.PIPE, stderr=sp.PIPE)
stdout, stderr = p.communicate()
temp1 = stderr.decode("utf-8")
temp2 = temp1.split("\n")
Rtemp, Ztemp, f1temp, f2temp, D2temp, Omegatemp, Utemp, Vtemp, veltemp = [] , [], [], [], [], [], [], [], []
for n1 in range(len(temp2)):
temp3 = temp2[n1].split(" ")
if temp3 == ['']:
pass
else:
Ztemp.append(float(temp3[0]))
Rtemp.append(float(temp3[1]))
f1temp.append(float(temp3[2]))
f2temp.append(float(temp3[3]))
D2temp.append(float(temp3[4]))
Omegatemp.append(float(temp3[5]))
Utemp.append(float(temp3[6]))
Vtemp.append(float(temp3[7]))
veltemp.append(float(temp3[8]))
R = np.asarray(Rtemp)
Z = np.asarray(Ztemp)
f1 = np.asarray(f1temp)
f2 = np.asarray(f2temp)
D2 = np.asarray(D2temp)
Omega = np.asarray(Omegatemp)
U = np.asarray(Utemp)
V = np.asarray(Vtemp)
vel = np.asarray(veltemp)
nr = int(len(R)/nz)
print("nr is %d" % nr)
R.resize((nz, nr))
Z.resize((nz, nr))
f1.resize((nz, nr))
f2.resize((nz, nr))
D2.resize((nz, nr))
Omega.resize((nz, nr))
U.resize((nz, nr))
V.resize((nz, nr))
vel.resize((nz, nr))
print('Got Field values')
return R, Z, f1, f2, D2, Omega, U, V, vel
# ------------------------------------------------------------------------------
nGFS = 2500
We = float(sys.argv[1])
Ohd= sys.argv[2]
Ohf = sys.argv[3]
hf = float(sys.argv[4])
print("Input Parameters: We = %f, Ohd = %s, Ohf = %s" % (We, Ohd, Ohf))
lw = 3
folder = 'TracerD2_InsideDrop' # output folder
if not os.path.isdir(folder):
os.makedirs(folder)
rmin, rmax, zmin, zmax = [-2.0, 2.0, -hf, 3.0-hf]
LEVEL = 10
nz = 2**(LEVEL)
for ti in range(nGFS):
t = 0.01 * ti
place = "intermediate/snapshot-%5.4f" % t
name = "%s/%8.8d.png" %(folder, int(t*1e3))
if not os.path.exists(place):
print("File %s not found!" % place)
else:
if os.path.exists(name):
print("Image %s found!" % name)
else:
facets1 = gettingFacets(place, 1)
if (len(facets1)):
facets2 = gettingFacets(place, 2)
R, Z, f1, f2, D2, Omega, U, V, vel = gettingfield(place)
speed = np.sqrt(U**2 + V**2)
fig, ax = plt.subplots()
fig.set_size_inches(19.20, 10.80)
rc('axes', linewidth=2)
# plt.xticks(fontsize=30)
# plt.yticks(fontsize=30)
ax.plot([0, 0], [zmin, zmax],'-.',color='grey',linewidth=lw)
# ax.plot([rmin, rmax], [0, 0],'-',color='grey',linewidth=lw/2)
ax.plot([rmin, rmin], [zmin, zmax],'-',color='black',linewidth=lw)
ax.plot([rmin, rmax], [zmin, zmin],'-',color='black',linewidth=lw)
ax.plot([rmin, rmax], [zmax, zmax],'-',color='black',linewidth=lw)
ax.plot([rmax, rmax], [zmin, zmax],'-',color='black',linewidth=lw)
## omega
cntrl1 = ax.imshow(D2, interpolation='bilinear', cmap="hot_r", origin='lower', extent=[0, -rmax, zmin, zmax], vmax = 1, vmin = -3)
## V
maxs = speed.max()
if maxs > 0.:
MaxVel = np.sqrt(We)
cntrl2 = ax.imshow(speed, interpolation='bilinear', cmap="Blues", origin='lower', extent=[0, rmax, zmin, zmax], vmax = MaxVel, vmin = 0.)
line_segments1 = LineCollection(facets1, linewidths=4, colors='green', linestyle='solid')
ax.add_collection(line_segments1)
line_segments2 = LineCollection(facets2, linewidths=4, colors='green', linestyle='solid')
ax.add_collection(line_segments2)
ax.set_xlabel(r'$\mathcal{R}$', fontsize=40)
ax.set_ylabel(r'$\mathcal{Z}$', fontsize=40)
ax.set_aspect('equal')
ax.set_xlim(rmin, rmax)
ax.set_ylim(zmin, zmax)
ax.set_title(r'$t/t_{\gamma} = %3.2f$' % t, fontsize=30)
# l, b, w, h = ax.get_position().bounds
# cb1 = fig.add_axes([l-0.16, b, 0.03, h])
# c1 = plt.colorbar(cntrl1,cax=cb1,orientation='vertical')
# c1.ax.tick_params(labelsize=30)
# c1.set_label(r'$log_{10}\left(2\mathcal{O}h\|D_{ij}\|\right)$',fontsize=40)
# c1.ax.yaxis.set_major_formatter(StrMethodFormatter('{x:,.0f}'))
# cb2 = fig.add_axes([l+w+0.01, b, 0.03, h])
# c2 = plt.colorbar(cntrl2,cax=cb2,orientation='vertical')
# c2.ax.tick_params(labelsize=30)
# c2.set_label(r'$\|V_i\|$',fontsize=40)
# c2.ax.yaxis.set_major_formatter(StrMethodFormatter('{x:,.1f}'))
# plt.show()
ax.axis('off')
plt.savefig(name, bbox_inches="tight")
plt.close()
else:
print("Problem in the available file %s" % place)
print(("Done %d of %d" % (ti+1, nGFS)))
