source: trunk/anuga_work/development/gareth/tests/runup/runuplot.py @ 8308

Last change on this file since 8308 was 8308, checked in by davies, 12 years ago

Added extrapolate_velocity_second_order option in config.py
and getting various tests working in gareth/tests

File size: 6.7 KB
Line 
1"""View results of runup.py
2"""
3#---------------
4# Import Modules
5#---------------
6import anuga
7import struct
8import numpy
9import scipy
10import pylab
11
12from Scientific.IO.NetCDF import NetCDFFile
13
14#--------------
15# Get variables
16#--------------
17fid = NetCDFFile('runup_v2.sww')
18x = fid.variables['x']
19x = x.getValue()
20y = fid.variables['y']
21y = y.getValue()
22elev = fid.variables['elevation']
23elev = elev.getValue()
24stage = fid.variables['stage']
25stage = stage.getValue()
26xmom = fid.variables['xmomentum']
27xmom = xmom.getValue()
28ymom = fid.variables['ymomentum']
29ymom = ymom.getValue()
30#xvel2 = fid.variables['centroid_xvelocity']
31#xvel2 = xvel2.getValue()
32#yvel2 = fid.variables['centroid_yvelocity']
33#yvel2 = yvel2.getValue()
34vols=fid.variables['volumes']
35vols=vols.getValue()
36
37# Calculate centroids
38x_cent=vols[:,0]*0.0
39y_cent=vols[:,0]*0.0
40
41for i in range(0,len(x_cent)):
42    x_cent[i]=(x[vols[i,0]]+x[vols[i,1]]+x[vols[i,2]])/3.0
43    y_cent[i]=(y[vols[i,0]]+y[vols[i,1]]+y[vols[i,2]])/3.0
44
45# Read in centroid velocities
46#xvel_cent = numpy.arange(0.0,len(x_cent)*xmom.shape[0]*1.0).reshape(xmom.shape[0],len(x_cent))
47#yvel_cent = numpy.arange(0.0,len(x_cent)*xmom.shape[0]*1.0).reshape(xmom.shape[0],len(x_cent))
48#
49#xfile=open('xvel.out','rb')
50#floatsize=struct.calcsize('f')
51#for i in range(0,xmom.shape[0]):
52#    for j in range(0,len(x_cent)):
53#        data = xfile.read(floatsize)
54#        num = struct.unpack('f', data)
55#        #print num[0], i*len(x_cent)+j
56#        xvel_cent[i,j] = num[0]
57#
58#xfile.close()
59#
60#yfile=open('yvel.out','rb')
61#floatsize2=struct.calcsize('f')
62#for i in range(0,xmom.shape[0]):
63#    for j in range(0,len(x_cent)):
64#        data2 = yfile.read(floatsize2)
65#        num2 = struct.unpack('f', data2)
66#        #print num[0], i*len(x_cent)+j
67#        yvel_cent[i,j] = num2[0] 
68#
69#yfile.close()
70
71#vel_cent=(xvel_cent**2+yvel_cent**2)**(0.5)
72## Read in centroid velocities from file as a long string
73##xvel2=[]
74##for line in file('xvel.txt'):
75##    line = line.rstrip('\n')
76##    xvel2.append(line)
77## Coerce values to array
78##xvel_cent = numpy.arange(0.0,len(x_cent)*xmom.shape[0]*1.0).reshape(xmom.shape[0],len(x_cent))
79##for i in range(0,xmom.shape[0]):
80##    for j in range(0,len(x_cent)):
81##        xvel_cent[i,j] = eval(xvel2[i*len(x_cent)+j])
82##        #print i, j, xvel_cent[i,j], xvel2[i*len(x_cent)+j], eval(xvel2[i*len(x_cent)+j])
83##
84##yvel2=[]
85##for line in file('yvel.txt'):
86#    line = line.rstrip('\n')
87#    yvel2.append(line)
88#
89#yvel_cent = numpy.arange(0.0,len(x_cent)*xmom.shape[0]*1.0).reshape(xmom.shape[0],len(x_cent))
90#
91#for i in range(0,xmom.shape[0]):
92#    for j in range(0,len(x_cent)):
93#        yvel_cent[i,j] = eval(yvel2[i*len(x_cent)+j])
94#
95
96#for i in vols.
97
98#-------------------
99# Calculate Velocity
100#-------------------
101xvel=xmom*0
102yvel=ymom*0
103for i in range(xmom.shape[0]):
104    xvel[i,:]=xmom[i,:]/(stage[i,:]-elev+1.0e-06)*(stage[i,:]-elev>1.0e-03)
105    yvel[i,:]=ymom[i,:]/(stage[i,:]-elev+1.0e-06)*(stage[i,:]-elev>1.0e-03)
106
107vel=(xvel**2+yvel**2)**0.5
108#------------------
109# Select line
110#------------------
111v=(y==0)
112
113#--------------------
114# Make plot animation
115#--------------------
116pylab.close() #If the plot is open, there will be problems
117pylab.ion()
118
119if True:
120    line, = pylab.plot( (x[v].min(),x[v].max()) ,(xvel[:,v].min(),xvel[:,v].max() ) )
121    for i in range(xmom.shape[0]):
122        line.set_xdata(x[v])
123        line.set_ydata(xvel[i,v])
124        pylab.draw()
125        pylab.plot( (0,1),(0,0), 'r' )
126        pylab.title(str(i)+'/200') # : velocity does not converge to zero' )
127        pylab.xlabel('x')
128        pylab.ylabel('Velocity (m/s)')
129
130    pylab.savefig('runup_x_velocities.png')
131
132#pylab.clf()
133#pylab.close()
134
135#------------------------------------------------
136# Maximum y velocities -- occurs in output step 3
137#------------------------------------------------
138#print yvel[3,:].max(), yvel[3,:].min()
139#highx=yvel[3,:].argmax()
140#v=(x==x[highx])
141#pylab.plot(yvel[3,v])
142#pylab.title('y-velocity is not always zero at the boundaries, e.g. x='+str(x[highx])+' , t=0.3s')
143#pylab.xlabel('y')
144#pylab.ylabel('Velocity (m/s)')
145#pylab.savefig('runup_y_velocities.png')
146
147pylab.clf()
148pylab.plot(x[y==0.5],stage[5,y==0.5])
149pylab.plot(x[y==0.5],stage[5,y==0.5],'o')
150pylab.plot(x[y==0.5],elev[y==0.5])
151pylab.xlabel('x (m)')
152pylab.ylabel('z (m)')
153pylab.title('Free surface and bed at y==0.5, time = 1.0 second')
154pylab.savefig('elev_1s_v2.png')
155
156pylab.clf()
157pylab.plot(x[y==0.5],xvel[5,y==0.5])
158pylab.plot(x[y==0.5],xvel[5,y==0.5],'o')
159pylab.title('Velocity at y==0.500, time = 1.0 second')
160pylab.xlabel('x (m)')
161pylab.ylabel('Velocity (m/s)')
162pylab.savefig('vel1d_1s_v2.png')
163
164#pylab.clf()
165#pylab.plot(x_cent[y_cent==y_cent[80]],xvel_cent[15,y_cent==y_cent[80]])
166#pylab.plot(x_cent[y_cent==y_cent[80]],xvel_cent[15,y_cent==y_cent[80]],'o')
167#pylab.title('Velocity at y==0.500, time = 3.0 second')
168#pylab.xlabel('x (m)')
169#pylab.ylabel('Velocity (m/s)')
170#pylab.savefig('vel1d_3s_v2.png')
171#-------------------------------------
172# Final velocities plot
173#-------------------------------------
174#pylab.clf()
175#pylab.quiver(x,y,xvel[200,:],yvel[200,:])
176#pylab.xlabel('x')
177#pylab.ylabel('y')
178#pylab.title('The maximum speed is '+ str(vel[200,:].max()) + ' m/s')
179#pylab.savefig('final_vel_field.png')
180#print vel[200,:].max()
181
182pylab.clf()
183pylab.scatter(x,y,c=elev,edgecolors='none', s=25)
184pylab.colorbar()
185#pylab.quiver(x_cent,y_cent,xvel_cent[15,:],yvel_cent[15,:])
186#pylab.title('The maximum speed is '+ str(vel_cent[15,:].max()) + ' m/s at time 3.0s')
187pylab.quiver(x,y,xvel[5,:],yvel[5,:])
188pylab.title('The maximum speed is '+ str(vel[5,:].max()) + ' m/s at time 1.0s')
189pylab.savefig('vel_1s_v2.png')
190
191pylab.clf()
192pylab.scatter(x,y,c=elev,edgecolors='none', s=25)
193pylab.colorbar()
194#pylab.quiver(x_cent,y_cent,xvel_cent[150,:],yvel_cent[150,:])
195#pylab.title('The maximum speed is '+ str(vel_cent[150,:].max()) + ' m/s at time 30.0s')
196pylab.quiver(x,y,xvel[150,:],yvel[150,:])
197pylab.title('The maximum speed is '+ str(vel[150,:].max()) + ' m/s at time 30.0s')
198pylab.savefig('vel_30s_v2.png')
199
200
201
202pylab.clf()
203pylab.plot(x[y==0.5],stage[150,y==0.5])
204pylab.plot(x[y==0.5],stage[150,y==0.5],'o')
205pylab.plot(x[y==0.5],elev[y==0.5])
206pylab.xlabel('x (m)')
207pylab.ylabel('z (m)')
208pylab.title('Free surface and bed at y==0.5, time = 30.0 second')
209pylab.savefig('elev_30s_v2.png')
210pylab.clf()
211pylab.plot(x[y==0.5],xvel[150,y==0.5])
212pylab.plot(x[y==0.5],xvel[150,y==0.5],'o')
213pylab.title('Velocity at y==0.500, time = 30.0 second')
214pylab.xlabel('x (m)')
215pylab.ylabel('Velocity (m/s)')
216pylab.savefig('vel1d_30s_v2.png')
217#pylab.clf()
218#pylab.plot(x_cent[y_cent==y_cent[80]],xvel_cent[150,y_cent==y_cent[80]])
219#pylab.plot(x_cent[y_cent==y_cent[80]],xvel_cent[150,y_cent==y_cent[80]],'o')
220#pylab.title('Velocity at y==0.500, time = 30.0 second')
221#pylab.xlabel('x (m)')
222#pylab.ylabel('Velocity (m/s)')
223#pylab.savefig('vel1d_30s_v2.png')
Note: See TracBrowser for help on using the repository browser.