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

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/channel_floodplain

Files:

• channel_floodplain1.py (modified) (6 diffs)
• plotme.py (modified) (3 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')