source: trunk/anuga_work/development/gareth/tests/parabolic/parabolaplot.py @ 8308

"""View results of runup.py
"""
#---------------
# Import Modules
#---------------
import anuga
import struct
import numpy
import scipy
import pylab

from Scientific.IO.NetCDF import NetCDFFile

#--------------
# Get variables
#--------------
fid = NetCDFFile('parabola_v2.sww')
x = fid.variables['x']
x = x.getValue()
y = fid.variables['y']
y = y.getValue()
elev = fid.variables['elevation']
elev = elev.getValue()
stage = fid.variables['stage']
stage = stage.getValue()
xmom = fid.variables['xmomentum']
xmom = xmom.getValue()
ymom = fid.variables['ymomentum']
ymom = ymom.getValue()
time = fid.variables['time']
time = time.getValue()

#xvel2 = fid.variables['centroid_xvelocity']
#xvel2 = xvel2.getValue()
#yvel2 = fid.variables['centroid_yvelocity']
#yvel2 = yvel2.getValue()
vols=fid.variables['volumes']
vols=vols.getValue()

# Calculate centroids
x_cent=vols[:,0]*0.0
y_cent=vols[:,0]*0.0

for i in range(0,len(x_cent)):
    x_cent[i]=(x[vols[i,0]]+x[vols[i,1]]+x[vols[i,2]])/3.0
    y_cent[i]=(y[vols[i,0]]+y[vols[i,1]]+y[vols[i,2]])/3.0

# Read in centroid velocities
xvel_cent = numpy.arange(0.0,len(x_cent)*xmom.shape[0]*1.0).reshape(xmom.shape[0],len(x_cent))
yvel_cent = numpy.arange(0.0,len(x_cent)*xmom.shape[0]*1.0).reshape(xmom.shape[0],len(x_cent))

xfile=open('xvel.out','rb')
floatsize=struct.calcsize('f')
for i in range(0,xmom.shape[0]):
    for j in range(0,len(x_cent)):
        data = xfile.read(floatsize)
        num = struct.unpack('f', data)
        #print num[0], i*len(x_cent)+j
        xvel_cent[i,j] = num[0] 

xfile.close()

yfile=open('yvel.out','rb')
floatsize2=struct.calcsize('f')
for i in range(0,xmom.shape[0]):
    for j in range(0,len(x_cent)):
        data2 = yfile.read(floatsize2)
        num2 = struct.unpack('f', data2)
        #print num[0], i*len(x_cent)+j
        yvel_cent[i,j] = num2[0]  

yfile.close()

vel_cent=(xvel_cent**2+yvel_cent**2)**(0.5)
# Read in centroid velocities from file as a long string
#xvel2=[]
#for line in file('xvel.txt'):
#    line = line.rstrip('\n')
#    xvel2.append(line)
# Coerce values to array
#xvel_cent = numpy.arange(0.0,len(x_cent)*xmom.shape[0]*1.0).reshape(xmom.shape[0],len(x_cent))
#for i in range(0,xmom.shape[0]):
#    for j in range(0,len(x_cent)):
#        xvel_cent[i,j] = eval(xvel2[i*len(x_cent)+j])
#        #print i, j, xvel_cent[i,j], xvel2[i*len(x_cent)+j], eval(xvel2[i*len(x_cent)+j])
#
#yvel2=[]
#for line in file('yvel.txt'):
#    line = line.rstrip('\n')
#    yvel2.append(line)
#
#yvel_cent = numpy.arange(0.0,len(x_cent)*xmom.shape[0]*1.0).reshape(xmom.shape[0],len(x_cent))
#
#for i in range(0,xmom.shape[0]):
#    for j in range(0,len(x_cent)):
#        yvel_cent[i,j] = eval(yvel2[i*len(x_cent)+j])
#

#for i in vols.

#-------------------
# Calculate Velocity
#-------------------
xvel=xmom*0
yvel=ymom*0
for i in range(xmom.shape[0]):
    xvel[i,:]=xmom[i,:]/(stage[i,:]-elev+1.0e-06)*(stage[i,:]-elev>1.0e-03)
    yvel[i,:]=ymom[i,:]/(stage[i,:]-elev+1.0e-06)*(stage[i,:]-elev>1.0e-03)

vel=(xvel**2+yvel**2)**0.5


#for i in range(xmim.shape[0]):
#    xvel[i,:]=xvel[i,:]*(stage[i,:]>elev+0.01)

#------------------
# Select line
#------------------
v=(y==0.0)
v2=(y_cent==y_cent[40])

#--------------------
# Make plot animation
#--------------------
pylab.close() #If the plot is open, there will be problems
pylab.ion()

if True:
    line, = pylab.plot( (x[v].min(),x[v].max()) ,(stage[:,v].min(),stage[:,v].max() ) )
    for i in range(700,xmom.shape[0]):
        line.set_xdata(x[v])
        line.set_ydata(stage[i,v])
        pylab.draw()
        #pylab.plot( (0,1),(0,0), 'r' )
        pylab.plot(x[v],elev[v]) 
        pylab.title(str(i)+'/200') # : velocity does not converge to zero' )
        pylab.xlabel('x')
        pylab.ylabel('stage (m/s)')

if True:
    pylab.clf()
 
    line, = pylab.plot( (x[v].min(),x[v].max()) ,(xvel[:,v].min(),xvel[:,v].max() ), 'ro' )
    for i in range(700,xmom.shape[0]):
        line.set_xdata(x[v])
        line.set_ydata(xvel[i,v])
        pylab.draw()
        pylab.plot( (x[v].min(),x[v].max()),(0,0), 'ro' )
        #pylab.plot(x[v],elev[v]) 
        pylab.title(str(i)+'/200') # : velocity does not converge to zero' )
        pylab.xlabel('x')
        pylab.ylabel('xvel (m/s)')

    #line, = pylab.plot( (x_cent[v2].min(),x_cent[v2].max()) ,(xvel_cent[:,v2].min(),xvel_cent[:,v2].max() ) )
    #for i in range(xmom.shape[0]):
    #    line.set_xdata(x_cent[v2])
    #    line.set_ydata(xvel_cent[i,v2])
    #    pylab.draw()
    #    pylab.plot( (x_cent[v2].min(),x_cent[v2].max()),(0,0), 'r' )
    #    #pylab.plot(x[v],elev[v]) 
    #    pylab.title(str(i)+'/200') # : velocity does not converge to zero' )
    #    pylab.xlabel('x')
    #    pylab.ylabel('xvel (m/s)')


D0=4.0
L=10.0
A=2.0
g=9.8
#t=time[30]

omega=numpy.sqrt(2*D0*g)/L
T= 2*numpy.pi/omega
ppp= (abs(x-4.*L/2.)).argmin()
ppp2= (abs(x-4.*L/4.)).argmin()

w = D0 + 2*A*D0/(L**2)*numpy.cos(omega*time)*( (x[ppp]-4.*L/2.) -A/2.*numpy.cos(omega*time))
w2 = D0 + 2*A*D0/(L**2)*numpy.cos(omega*time)*( (x[ppp2]-4.*L/2.) -A/2.*numpy.cos(omega*time))
w2 = w2*(w2>elev[ppp2])+elev[ppp2]*(w2<=elev[ppp2])

pylab.clf()
pylab.plot(time,w)
pylab.plot(time,stage[:,ppp])
pylab.savefig('Stage_centre_v2.png')
#       pylab.savefig('runup_x_velocities.png')
pylab.clf()
pylab.plot(time,w2)
pylab.plot(time,stage[:,ppp2])
pylab.savefig('Stage_centre_v3.png')


pltind=numpy.argmin(abs(time-T*3))

pylab.clf()
pylab.plot(x[v],xvel[pltind,v])
pylab.title('Velocity near to time=3T')
pylab.savefig('Vel_3T_v2.png')

pltind=numpy.argmin(abs(time-T*3.25))

pylab.clf()
pylab.plot(x[v],xvel[pltind,v])
pylab.title('Velocity near to time=3.25T')
pylab.savefig('Vel_3_5T_v2.png')
#pylab.clf()
#pylab.close()

#------------------------------------------------
# Maximum y velocities -- occurs in output step 3
#------------------------------------------------
#print yvel[3,:].max(), yvel[3,:].min()
#highx=yvel[3,:].argmax()
#v=(x==x[highx])
#pylab.plot(yvel[3,v])
#pylab.title('y-velocity is not always zero at the boundaries, e.g. x='+str(x[highx])+' , t=0.3s')
#pylab.xlabel('y')
#pylab.ylabel('Velocity (m/s)')
#pylab.savefig('runup_y_velocities.png')

#pylab.clf()
#pylab.plot(x[y==0.5],stage[15,y==0.5])
#pylab.plot(x[y==0.5],stage[15,y==0.5],'o')
#pylab.plot(x[y==0.5],elev[y==0.5])
#pylab.xlabel('x (m)')
#pylab.ylabel('z (m)')
#pylab.title('Free surface and bed at y==0.5, time = 3.0 second')
#pylab.savefig('elev_3s_vdense.png')
#pylab.clf()
#pylab.plot(x_cent[y_cent==0.500],xvel_cent[15,y_cent==0.500])
#pylab.plot(x_cent[y_cent==0.500],xvel_cent[15,y_cent==0.500],'o')
#pylab.title('Velocity at y==0.500, time = 3.0 second')
#pylab.xlabel('x (m)')
#pylab.ylabel('Velocity (m/s)')
#pylab.savefig('vel1d_3s_vdense.png')
#-------------------------------------
# Final velocities plot
#-------------------------------------
#pylab.clf()
#pylab.quiver(x,y,xvel[200,:],yvel[200,:])
#pylab.xlabel('x')
#pylab.ylabel('y')
#pylab.title('The maximum speed is '+ str(vel[200,:].max()) + ' m/s')
#pylab.savefig('final_vel_field.png')
#print vel[200,:].max()

#pylab.clf()
#pylab.scatter(x,y,c=elev,edgecolors='none', s=25)
#pylab.colorbar()
#pylab.quiver(x_cent,y_cent,xvel_cent[15,:],yvel_cent[15,:])
#pylab.title('The maximum speed is '+ str(vel_cent[15,:].max()) + ' m/s at time 3.0s')
#pylab.quiver(x,y,xvel[15,:],yvel[15,:])
#pylab.title('The maximum speed is '+ str(vel[15,:].max()) + ' m/s at time 3.0s')
#pylab.savefig('vel_3s_vdense.png')

#pylab.clf()
#pylab.scatter(x,y,c=elev,edgecolors='none', s=25)
#pylab.colorbar()
#pylab.quiver(x_cent,y_cent,xvel_cent[150,:],yvel_cent[150,:])
#pylab.title('The maximum speed is '+ str(vel_cent[150,:].max()) + ' m/s at time 30.0s')
#pylab.quiver(x,y,xvel[150,:],yvel[150,:])
#pylab.title('The maximum speed is '+ str(vel[150,:].max()) + ' m/s at time 30.0s')
#pylab.savefig('vel_30s_vdense.png')



#pylab.clf()
#pylab.plot(x[y==0.5],stage[150,y==0.5])
#pylab.plot(x[y==0.5],stage[150,y==0.5],'o')
#pylab.plot(x[y==0.5],elev[y==0.5])
#pylab.xlabel('x (m)')
#pylab.ylabel('z (m)')
#pylab.title('Free surface and bed at y==0.5, time = 30.0 second')
#pylab.savefig('elev_30s_vdense.png')
#pylab.clf()
#pylab.plot(x_cent[y_cent==0.500],xvel_cent[150,y_cent==0.500])
#pylab.plot(x_cent[y_cent==0.500],xvel_cent[150,y_cent==0.500],'o')
#pylab.title('Velocity at y==0.500, time = 30.0 second')
#pylab.xlabel('x (m)')
#pylab.ylabel('Velocity (m/s)')
##pylab.savefig('vel1d_30s_vdense.png')