Changeset 8353 for trunk/anuga_work/development/gareth/tests/runup

Timestamp:

Mar 7, 2012, 9:49:06 PM (13 years ago)

Author:

davies

Message:

Updates to balanced_dev and associated tests

Location:

trunk/anuga_work/development/gareth/tests/runup

Files:

• runup.py (modified) (5 diffs)
• runuplot.py (modified) (2 diffs)

trunk/anuga_work/development/gareth/tests/runup/runup.py

@@ -4 +4 @@
 #Import Modules
 #--------
-from sys import path
 import anuga
 
 import numpy
-
-import struct
-#import scipy
-
-#from Numeric import *
 
 from math import sin, pi, exp
@@ -23 +17 @@
 #from swb_domain import *
 #path.append('/home/gareth/storage/anuga_clean/anuga_jan12/trunk/anuga_work/development/gareth/balanced_basic')
-from balanced_basic import *
+#from balanced_basic import *
+from balanced_dev import *
 #---------
 #Setup computational domain
@@ -33 +28 @@
 domain.set_datadir('.')                          # Use current folder
 domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
-domain.set_store_vertices_uniquely(True)
+#domain.set_store_vertices_uniquely(True)
 #------------------
 # Define topography
@@ -66 +61 @@
 # Associate boundary tags with boundary objects
 #----------------------------------------------
-domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom':Br})
+domain.set_boundary({'left': Br, 'right': Bw, 'top': Br, 'bottom':Br})
 
 #------------------------------
@@ -80 +75 @@
     print domain.timestepping_statistics()
 
-#    momx=domain.quantities['xmomentum'].centroid_values
-#    momy=domain.quantities['ymomentum'].centroid_values
-#    #mom2=domain.quantities['xmomentum'].vertex_values
-#    dep1=(domain.quantities['stage'].centroid_values-domain.quantities['elevation'].centroid_values+1.0e-06)
-#    #dep2=(domain.quantities['stage'].vertex_values-domain.quantities['elevation'].vertex_values+1.0e-06)
-#    velx=momx/dep1
-#    vely=momy/dep1
-#    #vel2=mom2/dep2
-#    #print vel1.max(), vel2.max()
-#    #print vel1.min(), vel2.min()
-#    #xwrite.write(velx)
-#    #print >> xwrite, str(velx)
-#    #print >> ywrite, str(vely)
-#    for i in velx:
-#        data=struct.pack('f',i)
-#        xwrite.write(data)
-#
-#    for j in vely:
-#        data2=struct.pack('f',j)
-#        ywrite.write(data2)
-#
-#    #print >> xwrite, \n
-#    #numpy.savetxt("xvel.txt",velx)
-#    #numpy.savetxt("yvel.txt",vely)
-#    #velx.tofile(xwrite," ")
-#    #vely.tofile(ywrite," ")
-#
-#xwrite.close()
-#ywrite.close()
 
 print 'Finished'

trunk/anuga_work/development/gareth/tests/runup/runuplot.py

@@ -5 +5 @@
 #---------------
 import anuga
-import struct
 import numpy
 import scipy
-import pylab
+from matplotlib import pyplot as pyplot
+import util # Routines to read in and work with ANUGA output
 
-from Scientific.IO.NetCDF import NetCDFFile
+p2=util.get_output('runup_v2.sww', minimum_allowed_height=1.0e-03)
+p=util.get_centroids(p2, velocity_extrapolation=True)
 
-#--------------
-# Get variables
-#--------------
-fid = NetCDFFile('runup_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()
-#xvel2 = fid.variables['centroid_xvelocity']
-#xvel2 = xvel2.getValue()
-#yvel2 = fid.variables['centroid_yvelocity']
-#yvel2 = yvel2.getValue()
-vols=fid.variables['volumes']
-vols=vols.getValue()
+#p=util.get_output('runup_v2.sww', minimum_allowed_height=1.0e-03)
 
-# 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
 #------------------
 # Select line
 #------------------
-v=(y==0)
+py_central=p.y[scipy.argmin(abs(p.y-0.5))]
+v=(p.y==p.y[py_central])
 
 #--------------------
 # Make plot animation
 #--------------------
-pylab.close() #If the plot is open, there will be problems
-pylab.ion()
+pyplot.close() #If the plot is open, there will be problems
+pyplot.ion()
 
 if True:
-    line, = pylab.plot( (x[v].min(),x[v].max()) ,(xvel[:,v].min(),xvel[:,v].max() ) )
-    for i in range(xmom.shape[0]):
-        line.set_xdata(x[v])
-        line.set_ydata(xvel[i,v])
-        pylab.draw()
-        pylab.plot( (0,1),(0,0), 'r' )
-        pylab.title(str(i)+'/200') # : velocity does not converge to zero' )
-        pylab.xlabel('x')
-        pylab.ylabel('Velocity (m/s)')
+    line, = pyplot.plot( (p.x[v].min(),p.x[v].max()) ,(p.xvel[:,v].min(),p.xvel[:,v].max() ) )
+    for i in range(p.xmom.shape[0]):
+        line.set_xdata(p.x[v])
+        line.set_ydata(p.xvel[i,v])
+        pyplot.draw()
+        pyplot.plot( (0,1),(0,0), 'r' )
+        pyplot.title(str(i)+'/200') # : velocity does not converge to zero' )
+        pyplot.xlabel('x')
+        pyplot.ylabel('Velocity (m/s)')
 
-    pylab.savefig('runup_x_velocities.png')
+    pyplot.savefig('runup_x_velocities.png')
 
-#pylab.clf()
-#pylab.close()
+#pyplot.clf()
+#pyplot.close()
 
 #------------------------------------------------
@@ -139 +49 @@
 #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')
+#pyplot.plot(yvel[3,v])
+#pyplot.title('y-velocity is not always zero at the boundaries, e.g. x='+str(x[highx])+' , t=0.3s')
+#pyplot.xlabel('y')
+#pyplot.ylabel('Velocity (m/s)')
+#pyplot.savefig('runup_y_velocities.png')
 
-pylab.clf()
-pylab.plot(x[y==0.5],stage[5,y==0.5])
-pylab.plot(x[y==0.5],stage[5,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 = 1.0 second')
-pylab.savefig('elev_1s_v2.png')
+pyplot.clf()
+pyplot.plot(p.x[v],p.stage[5,v])
+pyplot.plot(p.x[v],p.stage[5,v],'o')
+pyplot.plot(p.x[v],p.elev[v])
+pyplot.xlabel('x (m)')
+pyplot.ylabel('z (m)')
+pyplot.title('Free surface and bed at y==0.5, time = 1.0 second')
+pyplot.savefig('elev_1s_v2.png')
 
-pylab.clf()
-pylab.plot(x[y==0.5],xvel[5,y==0.5])
-pylab.plot(x[y==0.5],xvel[5,y==0.5],'o')
-pylab.title('Velocity at y==0.500, time = 1.0 second')
-pylab.xlabel('x (m)')
-pylab.ylabel('Velocity (m/s)')
-pylab.savefig('vel1d_1s_v2.png')
+pyplot.clf()
+pyplot.plot(p.x[v],p.xvel[5,v])
+pyplot.plot(p.x[v],p.xvel[5,v],'o')
+pyplot.title('Velocity at y==0.500, time = 1.0 second')
+pyplot.xlabel('x (m)')
+pyplot.ylabel('Velocity (m/s)')
+pyplot.savefig('vel1d_1s_v2.png')
 
-#pylab.clf()
-#pylab.plot(x_cent[y_cent==y_cent[80]],xvel_cent[15,y_cent==y_cent[80]])
-#pylab.plot(x_cent[y_cent==y_cent[80]],xvel_cent[15,y_cent==y_cent[80]],'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_v2.png')
+#pyplot.clf()
+#pyplot.plot(x_cent[y_cent==y_cent[80]],xvel_cent[15,y_cent==y_cent[80]])
+#pyplot.plot(x_cent[y_cent==y_cent[80]],xvel_cent[15,y_cent==y_cent[80]],'o')
+#pyplot.title('Velocity at y==0.500, time = 3.0 second')
+#pyplot.xlabel('x (m)')
+#pyplot.ylabel('Velocity (m/s)')
+#pyplot.savefig('vel1d_3s_v2.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')
+#pyplot.clf()
+#pyplot.quiver(x,y,xvel[200,:],yvel[200,:])
+#pyplot.xlabel('x')
+#pyplot.ylabel('y')
+#pyplot.title('The maximum speed is '+ str(vel[200,:].max()) + ' m/s')
+#pyplot.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[5,:],yvel[5,:])
-pylab.title('The maximum speed is '+ str(vel[5,:].max()) + ' m/s at time 1.0s')
-pylab.savefig('vel_1s_v2.png')
+pyplot.clf()
+pyplot.scatter(p.x,p.y,c=p.elev,edgecolors='none', s=25)
+pyplot.colorbar()
+#pyplot.quiver(x_cent,y_cent,xvel_cent[15,:],yvel_cent[15,:])
+#pyplot.title('The maximum speed is '+ str(vel_cent[15,:].max()) + ' m/s at time 3.0s')
+pyplot.quiver(p.x,p.y,p.xvel[5,:],p.yvel[5,:])
+pyplot.title('The maximum speed is '+ str(p.vel[5,:].max()) + ' m/s at time 1.0s')
+pyplot.savefig('vel_1s_v2.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_v2.png')
+pyplot.clf()
+pyplot.scatter(p.x,p.y,c=p.elev,edgecolors='none', s=25)
+pyplot.colorbar()
+#pyplot.quiver(x_cent,y_cent,xvel_cent[150,:],yvel_cent[150,:])
+#pyplot.title('The maximum speed is '+ str(vel_cent[150,:].max()) + ' m/s at time 30.0s')
+pyplot.quiver(p.x,p.y,p.xvel[150,:],p.yvel[150,:])
+pyplot.title('The maximum speed is '+ str(p.vel[150,:].max()) + ' m/s at time 30.0s')
+pyplot.savefig('vel_30s_v2.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_v2.png')
-pylab.clf()
-pylab.plot(x[y==0.5],xvel[150,y==0.5])
-pylab.plot(x[y==0.5],xvel[150,y==0.5],'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_v2.png')
-#pylab.clf()
-#pylab.plot(x_cent[y_cent==y_cent[80]],xvel_cent[150,y_cent==y_cent[80]])
-#pylab.plot(x_cent[y_cent==y_cent[80]],xvel_cent[150,y_cent==y_cent[80]],'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_v2.png')
+pyplot.clf()
+pyplot.plot(p.x[v],p.stage[150,v])
+pyplot.plot(p.x[v],p.stage[150,v],'o')
+pyplot.plot(p.x[v],p.elev[v])
+pyplot.xlabel('x (m)')
+pyplot.ylabel('z (m)')
+pyplot.title('Free surface and bed at y==0.5, time = 30.0 second')
+pyplot.savefig('elev_30s_v2.png')
+pyplot.clf()
+pyplot.plot(p.x[v],p.xvel[150,v])
+pyplot.plot(p.x[v],p.xvel[150,v],'o')
+pyplot.title('Velocity at y==0.500, time = 30.0 second')
+pyplot.xlabel('x (m)')
+pyplot.ylabel('Velocity (m/s)')
+pyplot.savefig('vel1d_30s_v2.png')
+#pyplot.clf()
+#pyplot.plot(x_cent[y_cent==y_cent[80]],xvel_cent[150,y_cent==y_cent[80]])
+#pyplot.plot(x_cent[y_cent==y_cent[80]],xvel_cent[150,y_cent==y_cent[80]],'o')
+#pyplot.title('Velocity at y==0.500, time = 30.0 second')
+#pyplot.xlabel('x (m)')
+#pyplot.ylabel('Velocity (m/s)')
+#pyplot.savefig('vel1d_30s_v2.png')