Changeset 7403

Timestamp:

Aug 24, 2009, 2:17:07 PM (15 years ago)

Author:

ole

Message:

Played with Ted's V channel

File:

• anuga_work/development/ted_rigby/test_flow_in_vee_basic.py (modified) (5 diffs)

anuga_work/development/ted_rigby/test_flow_in_vee_basic.py

@@ -36 +36 @@
 width      = 20.
 side_slope = 1.0
-dx = dy    = 1          # Resolution: of points (elev) grid on both axes
+dx = dy    = 2          # Resolution: of points (elev) grid on both axes
 ref_flow   = 5.0        # Inflow at head of vee
 slope      = 1.0/1000
@@ -48 +48 @@
 domain = Domain(points, vertices, boundary)   
 domain.set_name('test_flow_in_vee')              # Output name
-
-
-#------------------------------------------------------------------------------
-# Setup initial conditions
-#------------------------------------------------------------------------------
-
-def topography(x, y):
-    z=-x * slope + abs((y-10)* side_slope)
-    return z
-
-domain.set_quantity('elevation', topography)  # Use function for elevation
-domain.set_quantity('friction', mannings_n)   # Constant friction of conc surface   
-domain.set_quantity('stage',
-                    expression='elevation')   # Dry initial condition
+print domain.statistics()
 
 
@@ -72 +59 @@
 normal_depth  = ( ref_flow*mannings_n* ( (2.0*( (1.0+(1.0/side_slope)**2.0)**0.5) )**0.6667 ) / (side_slope*(slope**0.5))  )**0.375
 area          = (normal_depth**2.0) *side_slope
-wet_perimeter = 2* normal_depth*((1+(side_slope)**2.0)**0.5)
+wet_perimeter = 2*normal_depth*((1+(side_slope)**2.0)**0.5)
 hyd_radius    = area/wet_perimeter
 flow_velocity = ref_flow/area
-top_width     = 2* side_slope*normal_depth
+top_width     = 2*side_slope*normal_depth
 
 # Compute xmomentum to input into the ANUGA inflow boundary
 # assuming momentum is applied uniformly to inflow boundary by ANUGA
 
-xmomentum     = ref_flow/top_width  # this is all very unrealistic but should replicate ref_flow in ANUGA
+#xmomentum = ref_flow/top_width  # this is all very unrealistic but should replicate ref_flow in ANUGA
+xmomentum = ref_flow/width       # This only applies where channel is flat. 
+
+print 'ref_flow', ref_flow
+print 'xmomentum', xmomentum
+print 'normal_depth', normal_depth
+print 'width', width
+
+
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+
+def topography(x, y):
+    z = -x * slope + abs((y-10)*side_slope)
+    
+    # Let both ends of domain be square so that the applied xmomentum 
+    # represents flow correctly 
+    for i in range(len(x)):
+        if x[i] < 20 or x[i] > 380:
+            z[i] = -x[i] * slope
+            
+    return z
+    
+    
+def initial_stage(x, y):
+    """Set constant depth based on calculated normal flow
+    """
+
+    w = -x*slope + normal_depth
+    return w
+    
+    
+domain.set_quantity('elevation', topography)  # Use function for elevation
+domain.set_quantity('friction', mannings_n)   # Constant friction of conc surface   
+domain.set_quantity('stage', initial_stage) 
+domain.set_quantity('xmomentum', xmomentum) # Initial momentum commensurate with normal flow
+
+
 
 #------------------------------------------------------------------------------
@@ -89 +115 @@
 Br = Reflective_boundary(domain)                       # Solid reflective wall on edge of sides
 Bi = Dirichlet_boundary([normal_depth,xmomentum,0.0])  # inflow at top
-Bo = Transmissive_boundary(domain)                     # momentum as above bdry
+Bo = Dirichlet_boundary([-length*slope+normal_depth,xmomentum,0.0])  # Outflow at bottom
+
 domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})
 
@@ -98 +125 @@
 for t in domain.evolve(yieldstep=100.0, finaltime=finaltime):
     if verbose :
-     print domain.timestepping_statistics()
+        print domain.timestepping_statistics()
+     
+     
+    # Check flow across flowline at 200m at every yield step
+    q=domain.get_flow_through_cross_section([[200.0,0.0],[200.0,width]])
+    print '90 degree flowline: ANUGA Q = %f, Ref_flow = %f' %(q, ref_flow)