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

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

Files:

• bedslope_discretization_compare (added)
• bedslope_discretization_compare/bedslope.R (added)
• channel_floodplain/channel_floodplain1.py (modified) (6 diffs)
• channel_floodplain/plotme.py (modified) (3 diffs)
• dam_break (added)
• dam_break/dam_break.py (added)
• dam_break/plotme.py (added)
• parabolic/parabolaplot.py (modified) (2 diffs)
• parabolic/parabolic.py (modified) (3 diffs)
• runup/runup.py (modified) (5 diffs)
• runup/runuplot.py (modified) (2 diffs)
• runup_sinusoid/runup_sinusoid.py (modified) (5 diffs)
• runup_sinusoid/runup_sinusoidplot.py (modified) (1 diff)
• util.py (modified) (5 diffs)

trunk/anuga_work/development/gareth/tests/channel_floodplain/channel_floodplain1.py

@@ -10 +10 @@
 import numpy
 from anuga.structures.inlet_operator import Inlet_operator
-
+from anuga import *
+#from swb_domain import domain
 #from anuga import *
 #from balanced_basic import *
@@ -23 +24 @@
 floodplain_width = 14.0 # Model domain width
 floodplain_slope = 1./300.
-chan_initial_depth = 0.8 # Initial depth of water in the channel
+chan_initial_depth = 0.65 # Initial depth of water in the channel
 chan_bankfull_depth = 1.0 # Bankfull depth of the channel
 chan_width = 10.0 # Bankfull width of the channel
 bankwidth = 2. # Width of the bank regions -- note that these protrude into the channel
 man_n=0.03 # Manning's n
-l0 = 2.009 # Length scale associated with triangle side length in channel (min_triangle area = 0.5*l0^2)
+l0 = 1.000 # Length scale associated with triangle side length in channel (min_triangle area = 0.5*l0^2)
 
 assert chan_width < floodplain_width, \
@@ -81 +82 @@
 
 
-domain.set_name('channel_floodplain1_test') # Output name
-domain.extrapolate_velocity_second_order=True
+domain.set_name('channel_floodplain1_bal_dev_lowvisc') # Output name
+#domain.set_store_vertices_uniquely(True)
+#domain.use_edge_limiter=False
+#domain.extrapolate_velocity_second_order=False
 #------------------------------------------------------------------------------
 #
@@ -139 +142 @@
 
 # Define inlet operator 
+flow_in_yval=100.0
 if True:
-    line1 = [ [floodplain_width/2. - chan_width/2., 0.],\
-              [floodplain_width/2. + chan_width/2., 0.] \
+    line1 = [ [floodplain_width/2. - chan_width/2., flow_in_yval],\
+              [floodplain_width/2. + chan_width/2., flow_in_yval] \
               ]
     Qin = 0.5*(floodplain_slope*(chan_width*chan_initial_depth)**2.*man_n**(-2.)\
@@ -191 +195 @@
 Br = anuga.Reflective_boundary(domain) # Solid reflective wall
 Bt = anuga.Transmissive_boundary(domain) # Transmissive boundary
-#Bout_sub = anuga.Dirichlet_boundary( \
-#        [-floodplain_length*floodplain_slope - chan_bankfull_depth + \
-#        chan_initial_depth, 0., 0.]) #An outflow boundary for subcritical steady flow
+
+
+Bout_sub = anuga.Dirichlet_boundary( \
+        [-floodplain_length*floodplain_slope - chan_bankfull_depth + \
+        chan_initial_depth, 0., 0.]) #An outflow boundary for subcritical steady flow
 
 def outflow_stage_boundary(t):
@@ -226 +232 @@
 #------------------------------------------------------------------------------
 
-for t in domain.evolve(yieldstep=1.0, finaltime=800.0):
+for t in domain.evolve(yieldstep=2.0, finaltime=3200.0):
     print domain.timestepping_statistics()
+    xx=domain.quantities['ymomentum'].centroid_values
+    dd=(domain.quantities['stage'].centroid_values - domain.quantities['elevation'].centroid_values)
+    dd=dd*(dd>0.)
+
+    tmp = xx/(dd+1.0e-06)*(dd>0.0)
+    print tmp.max(), tmp.argmax(), tmp.min(),  tmp.argmin()
+
+    # Compute flow through cross-section -- check that the inflow boundary condition is doing its job
+    # This also provides another useful steady-state check
+    if( numpy.floor(t/100.) == t/100. ):
+        print '#### COMPUTING FLOW THROUGH CROSS-SECTIONS########'
+        s1 = domain.get_flow_through_cross_section([[0., floodplain_length-300.0], [floodplain_width, floodplain_length-300.0]])
+        s2 = domain.get_flow_through_cross_section([[0., floodplain_length-1.0], [floodplain_width, floodplain_length-1.0]])
+        
+        print 'Cross sectional flows: ', s1, s2 
+        print '$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'
+
+    #vv = domain.get_flow_through_cross_section
     #print domain.quantities['ymomentum'].get_integral()
     #print 'Qin = ', Qin

trunk/anuga_work/development/gareth/tests/channel_floodplain/plotme.py

@@ -1 +1 @@
 import util
-import pylab
+from matplotlib import pyplot as pyplot
+#import pylab
 
 # Time-index to plot outputs from
 index=800
 
-p = util.get_output('channel_floodplain1_test.sww')
-v = (p.x>5.0)*(p.x<7.0)
+#p2 = util.get_output('channel_floodplain1_bal_dev.sww', minimum_allowed_height=0.01)
+#p=p2
+#p=util.get_centroids(p2, velocity_extrapolation=True)
+
+p2 = util.get_output('channel_floodplain1_bal_dev.sww', 0.01)
+p=util.get_centroids(p2, velocity_extrapolation=True)
+
+
+#p2 = util.get_output('channel_floodplain1_test_edge_lim.sww', minimum_allowed_height=0.01)
+#p2 = util.get_output('channel_floodplain1_test.sww', minimum_allowed_height=0.01)
+#p=util.get_centroids(p2, velocity_extrapolation=False)
+
+#p2 = util.get_output('channel_floodplain1_standard.sww', minimum_allowed_height=0.01)
+#p2 = util.get_output('channel_floodplain1_balanced_basic.sww', minimum_allowed_height=0.01)
+#p2 = util.get_output('channel_floodplain1_dev.sww')
+#p=util.get_centroids(p2, velocity_extrapolation=True)
+#p=p2
+v = (p.x>6.0)*(p.x<8.0)
 #util.animate_1D(p.time, p.stage[:, v], p.y[v])
 
@@ -18 +35 @@
 ##########################################################################
 
-Qin=6.6339  # Inflow discharge
+#Qin=6.6339  # Inflow discharge
+#Qin=5.2
+Qin=4.6932
+#Qin=4.693286 # Inflow discharge
 slp=1./300. # Floodplain slope (= water slope for steady uniform flow)
 man_n=0.03  # Manning's n
@@ -53 +73 @@
 # Analytical solution has U*abs(U)*n^2 / D^(4./3.) = Sf = bed slope
 # Hence, the following two variables should be equal -- I have checked that x velocities are fairly small
-pylab.clf()
-pylab.figure(figsize=(12.,8.))
-pylab.plot(p.y[v], (V2**2)**0.5,'o', label='computed velocity')
-pylab.plot(p.y[v], V2*0+k*dc_analytical**(2./3.), 'o', label='Analytical velocity')
-pylab.plot(p.y[v], V1, 'o',label='computed depth')
-pylab.plot(p.y[v], V1*0. + dc_analytical, 'o', label='analytical depth')
-#pylab.plot( ( (1./300.)*V1**(4./3.)*(1./0.03)**2.)**0.5,'o', label='Analytical velocity based on computed depth')
-pylab.title('Mid channel velocities and depths, vs analytical velocities and depths')
+pyplot.clf()
+pyplot.figure(figsize=(12.,8.))
+pyplot.plot(p.y[v], (V2**2)**0.5,'o', label='computed velocity')
+pyplot.plot(p.y[v], V2*0+k*dc_analytical**(2./3.), 'o', label='Analytical velocity')
+pyplot.plot(p.y[v], V1, 'o',label='computed depth')
+pyplot.plot(p.y[v], V1*0. + dc_analytical, 'o', label='analytical depth')
+#pyplot.plot( ( (1./300.)*V1**(4./3.)*(1./0.03)**2.)**0.5,'o', label='Analytical velocity based on computed depth')
+pyplot.title('Mid channel velocities and depths, vs analytical velocities and depths')
 # But in my tests, they are not equal
-#pylab.legend( ('computed velocity', 'Analytical velocity', 'computed depth', 'analytical depth'))
-pylab.savefig('trapz_velocity_downstream_l0_eq_2t.png')
+pyplot.legend( ('computed velocity', 'Analytical velocity', 'computed depth', 'analytical depth'), loc=4)
+pyplot.savefig('trapz_velocity_downstream_l0_eq_1_EL.png')
 
 # Plot velocity over the cross-section
-v1 = (p.y<500.0)&(p.y>470.0)
+v1 = (p.y<500.0)&(p.y>490.0)
 
-pylab.clf()
+pyplot.clf()
 analytical_stage = min(p.elev[v1]) + dc_analytical
 analytic_vel = ( (1./300.)*(analytical_stage-p.elev[v1])**(4./3.)*(1./0.03)**2.)**0.5
 analytic_vel = analytic_vel*(analytical_stage>p.elev[v1])
-pylab.figure(figsize=(12.,8.))
-pylab.plot(p.x[v1], p.yvel[index,v1],'o', label='computed velocity (m/s)')
-pylab.plot(p.x[v1], analytic_vel,'o', label='analytical velocity (m/s)')
-pylab.plot(p.x[v1],p.elev[v1],'o', label='bed elevation (m)')
-pylab.plot(p.x[v1],p.stage[index,v1],'o', label='computed stage (m)')
-pylab.plot(p.x[v1],p.stage[index,v1]*0. + analytical_stage,'o', label='analytical stage (m)')
+pyplot.figure(figsize=(12.,8.))
+pyplot.plot(p.x[v1], p.yvel[index,v1],'o', label='computed velocity (m/s)')
+pyplot.plot(p.x[v1], analytic_vel,'o', label='analytical velocity (m/s)')
+pyplot.plot(p.x[v1],p.elev[v1],'o', label='bed elevation (m)')
+pyplot.plot(p.x[v1],p.stage[index,v1],'o', label='computed stage (m)')
+pyplot.plot(p.x[v1],p.stage[index,v1]*0. + analytical_stage,'o', label='analytical stage (m)')
 
-#pylab.legend( ('computed velocity (m/s)', 'analytical velocity (m/s)', 'bed elevation (m)', 'computed stage (m)', 'analytical_stage (m)') )
-pylab.title('Velocity (analytical and numerical) and Stage:' + '\n' +'Central channel regions (470 to 500m)' +'\n')
-pylab.savefig('trapz_velocity_cross_channel_l0_eq_2t.png') 
-
+pyplot.legend( ('computed velocity (m/s)', 'analytical velocity (m/s)', 'bed elevation (m)', 'computed stage (m)', 'analytical_stage (m)') ,loc=10)
+pyplot.title('Velocity (analytical and numerical) and Stage:' + '\n' +'Central channel regions (470 to 500m)' +'\n')
+pyplot.savefig('trapz_velocity_cross_channel_l0_eq_1_EL.png') 
 
 
 # Plot velocity over the cross-section
-v1 = (p.y<800.0)&(p.y>770.0)
+v1 = (p.y<800.0)&(p.y>790.0)
 
-pylab.clf()
+pyplot.clf()
 analytical_stage = min(p.elev[v1]) + dc_analytical
 analytic_vel = ( (1./300.)*(analytical_stage-p.elev[v1])**(4./3.)*(1./0.03)**2.)**0.5
 analytic_vel = analytic_vel*(analytical_stage>p.elev[v1])
-pylab.figure(figsize=(12.,8.))
-pylab.plot(p.x[v1], p.yvel[index,v1],'o', label='computed velocity (m/s)')
-pylab.plot(p.x[v1], analytic_vel,'o', label='analytical velocity (m/s)')
-pylab.plot(p.x[v1],p.elev[v1],'o', label='bed elevation (m)')
-pylab.plot(p.x[v1],p.stage[index,v1],'o', label='computed stage (m)')
-pylab.plot(p.x[v1],p.stage[index,v1]*0. + analytical_stage,'o', label='analytical stage (m)')
+pyplot.figure(figsize=(12.,8.))
+pyplot.plot(p.x[v1], p.yvel[index,v1],'o', label='computed velocity (m/s)')
+pyplot.plot(p.x[v1], analytic_vel,'o', label='analytical velocity (m/s)')
+pyplot.plot(p.x[v1],p.elev[v1],'o', label='bed elevation (m)')
+pyplot.plot(p.x[v1],p.stage[index,v1],'o', label='computed stage (m)')
+pyplot.plot(p.x[v1],p.stage[index,v1]*0. + analytical_stage,'o', label='analytical stage (m)')
 
-#pylab.legend( ('computed velocity (m/s)', 'analytical velocity (m/s)', 'bed elevation (m)', 'computed stage (m)', 'analytical_stage (m)') )
-pylab.title('Velocity (analytical and numerical) and Stage:' + '\n' +'Central channel regions (770 to 800m)' +'\n')
-pylab.savefig('trapz_velocity_cross_channel_l0_eq_2tb.png') 
-
-
+pyplot.legend( ('computed velocity (m/s)', 'analytical velocity (m/s)', 'bed elevation (m)', 'computed stage (m)', 'analytical_stage (m)') ,loc=10)
+pyplot.title('Velocity (analytical and numerical) and Stage:' + '\n' +'Central channel regions (470 to 500m)' +'\n')
+pyplot.savefig('trapz_velocity_cross_channel_l0_eq_1b_EL.png') 
 
 # Plot velocity over the cross-section
-v1 = (p.y<900.0)&(p.y>870.0)
+v1 = (p.y<900.0)&(p.y>890.0)
 
-pylab.clf()
+pyplot.clf()
 analytical_stage = min(p.elev[v1]) + dc_analytical
 analytic_vel = ( (1./300.)*(analytical_stage-p.elev[v1])**(4./3.)*(1./0.03)**2.)**0.5
 analytic_vel = analytic_vel*(analytical_stage>p.elev[v1])
-pylab.figure(figsize=(12.,8.))
-pylab.plot(p.x[v1], p.yvel[index,v1],'o', label='computed velocity (m/s)')
-pylab.plot(p.x[v1], analytic_vel,'o', label='analytical velocity (m/s)')
-pylab.plot(p.x[v1],p.elev[v1],'o', label='bed elevation (m)')
-pylab.plot(p.x[v1],p.stage[index,v1],'o', label='computed stage (m)')
-pylab.plot(p.x[v1],p.stage[index,v1]*0. + analytical_stage,'o', label='analytical stage (m)')
+pyplot.figure(figsize=(12.,8.))
+pyplot.plot(p.x[v1], p.yvel[index,v1],'o', label='computed velocity (m/s)')
+pyplot.plot(p.x[v1], analytic_vel,'o', label='analytical velocity (m/s)')
+pyplot.plot(p.x[v1],p.elev[v1],'o', label='bed elevation (m)')
+pyplot.plot(p.x[v1],p.stage[index,v1],'o', label='computed stage (m)')
+pyplot.plot(p.x[v1],p.stage[index,v1]*0. + analytical_stage,'o', label='analytical stage (m)')
 
-#pylab.legend( ('computed velocity (m/s)', 'analytical velocity (m/s)', 'bed elevation (m)', 'computed stage (m)', 'analytical_stage (m)') )
-pylab.title('Velocity (analytical and numerical) and Stage:' + '\n' +'Central channel regions (870 to 900m)' +'\n')
-pylab.savefig('trapz_velocity_cross_channel_l0_eq_2tc.png') 
+pyplot.legend( ('computed velocity (m/s)', 'analytical velocity (m/s)', 'bed elevation (m)', 'computed stage (m)', 'analytical_stage (m)') , loc=10)
+pyplot.title('Velocity (analytical and numerical) and Stage:' + '\n' +'Central channel regions (870 to 900m)' +'\n')
+pyplot.savefig('trapz_velocity_cross_channel_l0_eq_1c_EL.png') 

trunk/anuga_work/development/gareth/tests/parabolic/parabolaplot.py

@@ -5 +5 @@
 #---------------
 import anuga
-import struct
+#import struct
 import numpy
-import scipy
-import pylab
+#import scipy
+from matplotlib import pyplot as pyplot
+import util
 
-from Scientific.IO.NetCDF import NetCDFFile
+p=util.get_output('parabola_v2.sww', 0.01)
+p2=util.get_centroids(p, velocity_extrapolation=True)
 
-#--------------
-# 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])
+# Define some 'lines' along which to plot
+v=(p2.y==p2.y[0])
 
 #--------------------
 # 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()) ,(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 False:
+    line, = pyplot.plot( (p2.x[v].min(),p2.x[v].max()) ,(p2.stage[:,v].min(),p2.stage[:,v].max() ) )
+    for i in range(700,p2.xmom.shape[0]):
+        line.set_xdata(p2.x[v])
+        line.set_ydata(p2.stage[i,v])
+        pyplot.draw()
+        #pyplot.plot( (0,1),(0,0), 'r' )
+        pyplot.plot(p2.x[v],p2.elev[v]) 
+        pyplot.title(str(i)+'/200') # : velocity does not converge to zero' )
+        pyplot.xlabel('x')
+        pyplot.ylabel('stage (m/s)')
 
-if True:
-    pylab.clf()
+if False:
+    pyplot.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)')
-
+    line, = pyplot.plot( (p2.x[v].min(),p2.x[v].max()) ,(p2.xvel[:,v].min(),p2.xvel[:,v].max() ), 'ro' )
+    for i in range(700,p2.xmom.shape[0]):
+        line.set_xdata(p2.x[v])
+        line.set_ydata(p2.xvel[i,v])
+        pyplot.draw()
+        pyplot.plot( (p2.x[v].min(),p2.x[v].max()),(0,0), 'ro' )
+        #pyplot.plot(x[v],elev[v]) 
+        pyplot.title(str(i)+'/200') # : velocity does not converge to zero' )
+        pyplot.xlabel('x')
+        pyplot.ylabel('xvel (m/s)')
 
 D0=4.0
@@ -172 +57 @@
 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()
+ppp= (abs(p2.x-4.*L/2.)).argmin()
+ppp2= (abs(p2.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])
+w = D0 + 2*A*D0/(L**2)*numpy.cos(omega*p2.time)*( (p2.x[ppp]-4.*L/2.) -A/2.*numpy.cos(omega*p2.time))
+w2 = D0 + 2*A*D0/(L**2)*numpy.cos(omega*p2.time)*( (p2.x[ppp2]-4.*L/2.) -A/2.*numpy.cos(omega*p2.time))
+w2 = w2*(w2>p2.elev[ppp2])+p2.elev[ppp2]*(w2<=p2.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')
+pyplot.clf()
+pyplot.plot(p2.time,w, color='blue')
+pyplot.plot(p2.time,p2.stage[:,ppp], color='green')
+pyplot.savefig('Stage_centre_v2.png')
+#       pyplot.savefig('runup_x_velocities.png')
+pyplot.clf()
+pyplot.plot(p2.time,w2, color='blue')
+pyplot.plot(p2.time,p2.stage[:,ppp2], color='green')
+pyplot.savefig('Stage_centre_v3.png')
 
 
-pltind=numpy.argmin(abs(time-T*3))
+pltind=numpy.argmin(abs(p2.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')
+pyplot.clf()
+pyplot.plot(p2.x[v],p2.xvel[pltind,v])
+pyplot.title('Velocity near to time=3T')
+pyplot.savefig('Vel_3T_v2.png')
 
-pltind=numpy.argmin(abs(time-T*3.25))
+pltind=numpy.argmin(abs(p2.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')
+pyplot.clf()
+pyplot.plot(p2.x[v],p2.xvel[pltind,v])
+pyplot.title('Velocity near to time=3.25T')
+pyplot.savefig('Vel_3_5T_v2.png')
+#pyplot.clf()
+#pyplot.close()

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

@@ -5 +5 @@
 #--------
 import anuga
-
+#from anuga.shallow_water.shallow_water_domain import Domain as Domain
 import numpy
 
-import struct
-#import scipy
-
-#from Numeric import *
-
-from math import sin, pi, exp
+from balanced_dev import *
+#from balanced_basic import *
 #from anuga.shallow_water_balanced2.swb2_domain import Domain as Domain
-from anuga.shallow_water.shallow_water_domain import Domain as Domain
+#from anuga.shallow_water.shallow_water_domain import Domain as Domain
 #---------
 #Setup computational domain
@@ -25 +21 @@
 domain.set_datadir('.')                          # Use current folder
 domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
-#domain.set_minimum_allowed_height(0.01)
+domain.set_minimum_allowed_height(0.01)
 # Time stepping
 #domain.set_timestepping_method('euler') # Default
@@ -81 +77 @@
 #Evolve the system through time
 #------------------------------
-
-xwrite=open("xvel.out","wb")
-ywrite=open("yvel.out","wb")
-# Set print options to be compatible with file writing via the 'print' statement 
-numpy.set_printoptions(threshold=numpy.nan, linewidth=numpy.nan)
-
 for t in domain.evolve(yieldstep=0.1,finaltime=90.0):
     print domain.timestepping_statistics()
-    # Calculate velocity at centroids, and record
-    momx=domain.quantities['xmomentum'].centroid_values
-    momy=domain.quantities['ymomentum'].centroid_values
-    dep1=(domain.quantities['stage'].centroid_values-domain.quantities['elevation'].centroid_values+1.0e-06)
-    velx=momx/dep1
-    vely=momy/dep1
-    
-    for i in velx:
-        data=struct.pack('f',i)
-        xwrite.write(data)
 
-    for j in vely:
-        data2=struct.pack('f',j)
-        ywrite.write(data2)
-
-
-xwrite.close()
-ywrite.close()
 print 'Finished'

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')

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

@@ -8 +8 @@
 import numpy
 
-import struct
-#import scipy
-from sys import path
-#from Numeric import *
-
 from math import sin, pi, exp
 #from anuga.shallow_water.shallow_water_domain import Domain as Domain
@@ -18 +13 @@
 #path.append('/home/gareth/storage/anuga_clean/anuga_jan12/trunk/anuga_work/development/gareth/balanced_basic')
 #from swb2_domain import *
-from balanced_basic import *
+#from balanced_basic import *
+from balanced_dev import *
 #---------
 #Setup computational domain
@@ -28 +24 @@
 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)
-# Time stepping
-#domain.set_timestepping_method('euler') # Default
-#domain.set_timestepping_method('rk2') # 
-#domain.beta_w= 1. #0.2
-#domain.beta_uh= 1. #0.2
-#domain.beta_vh= 1. #0.2
-#domain.beta_w_dry= 0.2 #0.
-#domain.beta_uh_dry= 0.2 #0.
-#domain.beta_vh_dry= 0.2 #0.
-#domain.alpha=100.
+#domain.set_store_vertices_uniquely(True)
 
 #------------------
@@ -45 +31 @@
 
 def topography(x,y):
-        return -x/2.0 +0.05*numpy.sin((x+y)*50.0) #+0.1*(numpy.random.rand(len(x)) -0.5)       # Linear bed slope + small random perturbation 
-# Note that without smoothing, this will cause the elevation to be discontinuous, which at present the model cannot manage.
+        return -x/2.0 +0.05*numpy.sin((x+y)*50.0)  
 
 def stagefun(x,y):
@@ -88 +73 @@
 #------------------------------
 
-xwrite=open("xvel.out","wb")
-ywrite=open("yvel.out","wb")
-# Set print options to be compatible with file writing via the 'print' statement 
-numpy.set_printoptions(threshold=numpy.nan, linewidth=numpy.nan)
-
 for t in domain.evolve(yieldstep=0.1,finaltime=20.0):
     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
-    dep1=dep1*(dep1>1.0e-03) +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_sinusoid/runup_sinusoidplot.py

@@ -8 +8 @@
 import numpy
 
-p1=util.get_output('runup_sinusoid_v2.sww')
-p2=util.get_centroids(p1)
+p1=util.get_output('runup_sinusoid_v2.sww', 0.001)
+p2=util.get_centroids(p1, velocity_extrapolation=True)
 
 #------------------

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

@@ -4 +4 @@
 class get_output:
     """Read in data from an .sww file in a convenient form
-       e.g. p = util.get_output('channel3.sww')
+       e.g. 
+        p = util.get_output('channel3.sww', minimum_allowed_height=0.01)
+        
+       p then contains most relevant information as e.g., p.stage, p.elev, p.xmom, etc 
     """
-    def __init__(self, filename):
+    def __init__(self, filename, minimum_allowed_height=1.0e-03):
         self.x, self.y, self.time, self.vols, self.stage, \
                 self.elev, self.xmom, self.ymom, \
-                self.xvel, self.yvel, self.vel = \
-                read_output(filename)
+                self.xvel, self.yvel, self.vel, self.minimum_allowed_height = \
+                read_output(filename, minimum_allowed_height)
 
 
-def read_output(filename):
+def read_output(filename, minimum_allowed_height):
     # Input: The name of an .sww file to read data from,
     #                    e.g. read_sww('channel3.sww')
@@ -61 +64 @@
 
     for i in range(xmom.shape[0]):
-        xvel[i,:]=xmom[i,:]/(stage[i,:]-elev+1.0e-06)*(stage[i,:]> elev+2.0e-02)
-        yvel[i,:]=ymom[i,:]/(stage[i,:]-elev+1.0e-06)*(stage[i,:]> elev+2.0e-02)
+        xvel[i,:]=xmom[i,:]/(stage[i,:]-elev+1.0e-06)*(stage[i,:]> elev+minimum_allowed_height)
+        yvel[i,:]=ymom[i,:]/(stage[i,:]-elev+1.0e-06)*(stage[i,:]> elev+minimum_allowed_height)
 
     vel = (xvel**2+yvel**2)**0.5
 
-    return x, y, time, vols, stage, elev, xmom, ymom, xvel, yvel, vel
+    return x, y, time, vols, stage, elev, xmom, ymom, xvel, yvel, vel, minimum_allowed_height
 
 ##############
@@ -102 +105 @@
 class get_centroids:
     """
-    Extract centroid values from p=util.get_output('my_sww.sww')
+    Extract centroid values from the output of get_output.  
+    e.g.
+        p = util.get_output('my_sww.sww', minimum_allowed_height=0.01) # vertex values
+        pc=util.get_centroids(p, velocity_extrapolation=True) # centroid values
     """
-    def __init__(self,p):
-        self.x, self.y, self.stage, self.xmom,\
+    def __init__(self,p, velocity_extrapolation=False):
+        self.time, self.x, self.y, self.stage, self.xmom,\
              self.ymom, self.elev, self.xvel, \
-             self.yvel, self.vel= get_centroid_values(p)
+             self.yvel, self.vel= \
+             get_centroid_values(p, velocity_extrapolation)
+                                 
 
-def get_centroid_values(p):
+def get_centroid_values(p, velocity_extrapolation):
     # Input: p is the result of e.g. p=util.get_output('mysww.sww'). See the get_output class defined above
     # Output: Values of x, y, Stage, xmom, ymom, elev, xvel, yvel, vel at centroids
@@ -120 +128 @@
     vols2=numpy.zeros(l, dtype='int')
 
+    # FIXME: 22/2/12/ - I think this loop is slow, should be able to do this
+    # another way
     for i in range(l):
         vols0[i]=p.vols[i][0]
@@ -130 +140 @@
 
     stage_cent=(p.stage[:,vols0]+p.stage[:,vols1]+p.stage[:,vols2])/3.0
-    xmom_cent=(p.xmom[:,vols0]+p.xmom[:,vols1]+p.xmom[:,vols2])/3.0
-    ymom_cent=(p.ymom[:,vols0]+p.ymom[:,vols1]+p.ymom[:,vols2])/3.0
     elev_cent=(p.elev[vols0]+p.elev[vols1]+p.elev[vols2])/3.0
 
-    # Now compute velocities
-    xvel_cent=stage_cent*0.0
-    yvel_cent=stage_cent*0.0
-    vel_cent=stage_cent*0.0
+    # Here, we need to treat differently the case of momentum extrapolation or
+    # velocity extrapolation
+    if velocity_extrapolation:
+        xvel_cent=(p.xvel[:,vols0]+p.xvel[:,vols1]+p.xvel[:,vols2])/3.0
+        yvel_cent=(p.yvel[:,vols0]+p.yvel[:,vols1]+p.yvel[:,vols2])/3.0
+        
+        # Now compute momenta
+        xmom_cent=stage_cent*0.0
+        ymom_cent=stage_cent*0.0
 
-    t=len(p.time)
+        t=len(p.time)
 
-    for i in range(t):
-        xvel_cent[i,:]=xmom_cent[i,:]/(stage_cent[i,:]-elev_cent+1.0e-06)*(stage_cent[i,:]>elev_cent+1.0e-03)
-        yvel_cent[i,:]=ymom_cent[i,:]/(stage_cent[i,:]-elev_cent+1.0e-06)*(stage_cent[i,:]>elev_cent+1.0e-03)
-    
+        for i in range(t):
+            xmom_cent[i,:]=xvel_cent[i,:]*(stage_cent[i,:]-elev_cent+1.0e-06)*\
+                                            (stage_cent[i,:]>elev_cent+p.minimum_allowed_height)
+            ymom_cent[i,:]=yvel_cent[i,:]*(stage_cent[i,:]-elev_cent+1.0e-06)*\
+                                            (stage_cent[i,:]>elev_cent+p.minimum_allowed_height)
+
+    else:
+        xmom_cent=(p.xmom[:,vols0]+p.xmom[:,vols1]+p.xmom[:,vols2])/3.0
+        ymom_cent=(p.ymom[:,vols0]+p.ymom[:,vols1]+p.ymom[:,vols2])/3.0
+
+        # Now compute velocities
+        xvel_cent=stage_cent*0.0
+        yvel_cent=stage_cent*0.0
+
+        t=len(p.time)
+
+        for i in range(t):
+            xvel_cent[i,:]=xmom_cent[i,:]/(stage_cent[i,:]-elev_cent+1.0e-06)*(stage_cent[i,:]>elev_cent+p.minimum_allowed_height)
+            yvel_cent[i,:]=ymom_cent[i,:]/(stage_cent[i,:]-elev_cent+1.0e-06)*(stage_cent[i,:]>elev_cent+p.minimum_allowed_height)
+   
+
+    # Compute velocity 
     vel_cent=(xvel_cent**2 + yvel_cent**2)**0.5
 
-    return x_cent, y_cent, stage_cent, xmom_cent,\
+    return p.time, x_cent, y_cent, stage_cent, xmom_cent,\
              ymom_cent, elev_cent, xvel_cent, yvel_cent, vel_cent