Plotting

In [ ]:

    

import numpy as np
import matplotlib.pyplot as plt
% matplotlib inline
from matplotlib import animation, cm
from IPython.display import HTML
from numpy.linalg import inv
from mpl_toolkits.mplot3d import Axes3D
import tables as tb


    

In [ ]:

    

# read input file 
input_file = open('input_file.txt', 'r')
inputs = input_file.readlines()

for line in inputs:
    name, *dat = line.split();
    
    if not (name == 'periodic' or name == 'outfile'):
        if len(dat) == 1:
            if name in ('nx', 'ny', 'nt', 'nlayers', 'dprint'):
                exec(name + ' = int(dat[0])')
            else:
                    
                exec(name + ' = float(dat[0])')
        else:
            exec(name + ' = [float(d) for d in dat]')
    elif name == 'outfile':
        outfile = dat[0]

            
dx = (xmax - xmin) / (nx-2)
dy = (ymax - ymin) / (ny-2)
dt = 0.1 * min(dx, dy)
input_file.close()


    

In [ ]:

    

# read data
outfile = '../../Documents/Work/swerve/out.h5'
#outfile = 'out.dat'
if (outfile[-2:] == 'h5'): #hdf5
    f = tb.open_file(outfile, 'r')
    table = f.root.SwerveOutput
    D_2d = np.swapaxes(table[:,:,:,:,0], 1, 3)
else: # assume some kind of csv
    nx      =150
    ny      =150
    nt      =600
    nlayers =2
    xmin   = 0.0
    xmax    =10.0
    ymin    =0.0
    ymax    =10.0
    data = np.loadtxt(outfile, delimiter=',')
    ts = data[:,0]
    xs = data[:,1]
    ys = data[:,2]
    ls = data[:,3]
    Ds = data[:,4]
    #Sxs = data[:,5]
    #Syx = data[:,6]
    #t = range(int((nt+1)/dprint))*dprint
    #print(len(Ds))
    
    D_2d = np.zeros((nt, nlayers, nx, ny))

    #D_2d = np.zeros((len(Ds), nlayers, nx, ny))
    #Sx_2d = np.zeros((nt, nlayers, nx, ny))
    #Sy_2d = np.zeros((nt, nlayers, nx, ny))

    for i in range(len(Ds)):
        #print(int(xs[i]*nx/xmax))
        #print(int(ts[i]), int(ls[i]), int(xs[i]), int(ys[i]))
        D_2d[int(ts[i]), int(ls[i]), int(xs[i]), int(ys[i])] = Ds[i]


    

In [ ]:

    

x = np.linspace(0, xmax, num=nx-4, endpoint=False)
y = np.linspace(0, ymax, num=ny-4, endpoint=False)


    

In [ ]:

    

X, Y = np.meshgrid(x,y)


    

In [ ]:

    

print(np.shape(X))
print(np.shape(Y))
print(np.shape(D_2d[0,1,2:-2,2:-2].T))


    

In [ ]:

    

fig = plt.figure(figsize=(12,10))
ax = fig.gca(projection='3d')

ax.set_xlim(0,10)
ax.set_ylim(0,10)
ax.set_zlim(0.7,1.4)
#ax.set_zlim(-2,5)

ax.plot_surface(X,Y,D_2d[1,1,2:-2,2:-2].T, rstride=1, cstride=2, lw=0, cmap=cm.viridis, antialiased=True)

ax.plot_wireframe(X,Y,D_2d[1,0,2:-2,2:-2].T, rstride=2, cstride=2, lw=0.1, cmap=cm.viridis, antialiased=True)

#plt.plot(x,h[0,1:-1,0],x,h[1,1:-1,0], lw=2)
plt.show()


    

In [ ]:

    

fig = plt.figure(figsize=(12,10))
ax = fig.gca(projection='3d')

surface_1 = ax.plot_surface(X,Y,D_2d[0,1,2:-2,2:-2].T, rstride=1, cstride=2, lw=0, cmap=cm.viridis, antialiased=True)
surface_2 = ax.plot_wireframe(X,Y,D_2d[0,0,2:-2,2:-2].T, rstride=2, cstride=2, lw=0.1, cmap=cm.viridis, antialiased=True)

def animate(i):
    ax.clear()
    ax.set_xlim(0,10)
    ax.set_ylim(0,10)
    ax.set_zlim(0.7,1.4)
    #ax.view_init(80)
    ax.plot_surface(X,Y,D_2d[i,1,2:-2,2:-2].T, rstride=1, cstride=2, lw=0, cmap=cm.viridis, antialiased=True)
    ax.plot_wireframe(X,Y,D_2d[i,0,2:-2,2:-2].T, rstride=2, cstride=2, lw=0.1, cmap=cm.viridis, antialiased=True)

anim = animation.FuncAnimation(fig, animate, frames=50, interval=60)


    

In [ ]:

    

HTML(anim.to_html5_video())


    

In [ ]:

    

to_save = False

if to_save:
    Writer = animation.writers['ffmpeg']
    writer = Writer(fps=5, metadata=dict(artist='Me'), bitrate=1800)
    anim.save('gr_cuda.mp4', writer=writer)


    

In [ ]:

    

if (outfile[-2:] == 'h5'): #hdf5
    f.close()