Changeset 8196

Timestamp:

Aug 12, 2011, 1:30:07 PM (13 years ago)

Author:

paul

Message:

Pipe now has two reconstruction methods available, both taken from channel

Location:

trunk/anuga_work/development/2010-projects/anuga_1d/pipe

Files:

• pipe_domain.py (modified) (4 diffs)
• pipe_domain_ext.c (modified) (2 diffs)
• test_pipe.py (modified) (3 diffs)

trunk/anuga_work/development/2010-projects/anuga_1d/pipe/pipe_domain.py

@@ -140 +140 @@
         #Call correct module function
         #(either from this module or C-extension)
-        distribute_to_vertices_and_edges_limit_a_d(self)
+        distribute_to_vertices_and_edges_limit_s_v_h(self)
+        #distribute_to_vertices_and_edges_limit_s_v(self)
         
 
@@ -166 +167 @@
 
 #-----------------------------------------------------------------------
-# Distribute to verticies with stage reconstructed and then extrapolated
+# Distribute to verticies with stage, velocity and pipe geometry
+# reconstructed and then extrapolated.
 #-----------------------------------------------------------------------
-def distribute_to_vertices_and_edges_limit_a_d(domain):
-    
-    #Remove very thin layers of water
-    #protect_against_infinitesimal_and_negative_heights(domain)  
-
-
-
+def distribute_to_vertices_and_edges_limit_s_v(domain):
     import sys
     from anuga_1d.config import epsilon, h0
 
-    
     N = domain.number_of_elements
 
@@ -203 +198 @@
     s_C   = state.centroid_values
 
-    if domain.setstageflag:
-        a_C[:,] = (w_C-z_C)*b_C
-        domain.setstageflag = False
-
-    # Paul: Same for top?
-    if domain.discontinousb:
-        width.extrapolate_second_order()
-        
-
-    h0 = 1.0e-12
-
-    # Paul: How does top (i.e. t_C) affect these?
-    h_C[:] = numpy.where( (a_C>h0) | (b_C>h0), a_C/(b_C+h0/b_C), 0.0 )
-    u_C[:] = numpy.where( (a_C>h0) | (b_C>h0), d_C/(a_C+h0/a_C), 0.0 )
-
-    a_C[:] = numpy.where( (a_C>h0) | (b_C>h0), a_C, 0.0 )
-    b_C[:] = numpy.where( (a_C>h0) | (b_C>h0), b_C, 0.0 )
-    t_C[:] = numpy.where( (a_C>h0) | (b_C>h0), t_C, 0.0 )
-    d_C[:] = numpy.where( (a_C>h0) | (b_C>h0), d_C, 0.0 )
-
-    w_C[:] = h_C + z_C
-    
-                
-    # Paul: same for top?
-    bed.extrapolate_second_order()
-    
-    for name in ['velocity','stage']:
+    # Calculate height, velocity and stage.
+    # Here we assume the conserved quantities and the channel geometry
+    # (i.e. bed, top and width) have been accurately computed in the previous
+    # timestep. 
+    h_C[:] = numpy.where(a_C > 0.0, a_C/(b_C + h0/b_C), 0.0)
+    u_C[:] = numpy.where(a_C > 0.0, d_C/(a_C + h0/a_C), 0.0)
+
+    w_C[:] = h_C + z_C    
+      
+    # Extrapolate velocity and stage as well as channel geometry.
+    for name in ['velocity', 'stage', 'elevation', 'width', 'top']:
         Q = domain.quantities[name]
         if domain.order == 1:
@@ -238 +217 @@
             raise 'Unknown order'
 
-
-    a_V  = area.vertex_values
+    # Stage, bed, width, top and velocity have been extrapolated.
+    # We may need to ensure that top is never 0.
     w_V  = stage.vertex_values
+    u_V  = velocity.vertex_values
     z_V  = bed.vertex_values
-    h_V  = height.vertex_values
-    u_V  = velocity.vertex_values
-    d_V  = discharge.vertex_values
     b_V  = width.vertex_values
     t_V  = top.vertex_values
+
+    # These quantites need to update vertex_values
+    a_V  = area.vertex_values
+    h_V  = height.vertex_values
+    d_V  = discharge.vertex_values
+
+    # State is true of false so should not be extrapolated
     s_V  = state.vertex_values
 
 
-    h_V[:] = w_V-z_V
-
-    w_V[:] = numpy.where(h_V > h0, w_V, z_V)
-    h_V[:] = numpy.where(h_V > h0, h_V, 0.0)
+    # Calculate height and fix up negatives. The main idea here is
+    # fix up the wet/dry interface.
+    h_V[:,:] = w_V - z_V
+
+    h_0 = numpy.where(h_V[:,0] < 0.0, 0.0, h_V[:,0])
+    h_1 = numpy.where(h_V[:,0] < 0.0, h_V[:,1]+h_V[:,0], h_V[:,1])
+
+    h_V[:,0] = h_0
+    h_V[:,1] = h_1
+
+
+    h_0 = numpy.where(h_V[:,1] < 0.0, h_V[:,1]+h_V[:,0], h_V[:,0])
+    h_1 = numpy.where(h_V[:,1] < 0.0, 0.0, h_V[:,1])
+
+    h_V[:,0] = h_0
+    h_V[:,1] = h_1
+
+
+    # Protect against negative and small heights. If we set h to zero
+    # we better do the same with velocity (i.e. no water, no velocity).
+    h_V[:,:] = numpy.where (h_V <= h0, 0.0, h_V)
+    u_V[:,:] = numpy.where (h_V <= h0, 0.0, u_V)
+
+
+    # Clean up conserved quantities
+    w_V[:] = z_V + h_V
+    a_V[:] = b_V * h_V
+    d_V[:] = u_V * a_V
+
+    
+    return
+
+#-----------------------------------------------------------------------
+# Distribute to verticies with stage, height and velocity reconstructed 
+# and then extrapolated.
+# In this method, we extrapolate the stage and height out to the vertices.
+# The bed, although given as initial data to the problem, is reconstructed
+# from the stage and height. This ensures consistency of the reconstructed
+# quantities (i.e. w = z + h) as well as protecting against negative
+# heights.
+#-----------------------------------------------------------------------
+def distribute_to_vertices_and_edges_limit_s_v_h(domain):
+    import sys
+    from anuga_1d.config import epsilon, h0
+
+    N = domain.number_of_elements
+
+    #Shortcuts
+    area      = domain.quantities['area']
+    discharge = domain.quantities['discharge']
+    bed       = domain.quantities['elevation']
+    height    = domain.quantities['height']
+    velocity  = domain.quantities['velocity']
+    width     = domain.quantities['width']
+    top     = domain.quantities['top']
+    stage     = domain.quantities['stage']
+    state     = domain.quantities['state']
+
+    #Arrays   
+    a_C   = area.centroid_values
+    d_C   = discharge.centroid_values
+    z_C   = bed.centroid_values
+    h_C   = height.centroid_values
+    u_C   = velocity.centroid_values
+    b_C   = width.centroid_values
+    t_C   = top.centroid_values
+    w_C   = stage.centroid_values
+    s_C   = state.centroid_values
+
+    # Construct h,u,w from the conserved quantities after protecting
+    # conserved quantities from becoming too small.
+    a_C[:] = numpy.where( (a_C>h0), a_C, 0.0 )
+    d_C[:] = numpy.where( (a_C>h0), d_C, 0.0 )
+    h_C[:] = numpy.where( (b_C>h0), a_C/(b_C + h0/b_C), 0.0 )
+    u_C[:] = numpy.where( (a_C>h0), d_C/(a_C + h0/a_C), 0.0 )    
+
+    # Set the stage
+    w_C[:] = h_C + z_C          
+
+    # Extrapolate "fundamental" quantities.
+    # All other quantities will be reconstructed from these.
+    for name in ['velocity', 'stage',  'height', 'width', 'top']:
+        Q = domain.quantities[name]
+        if domain.order == 1:
+            Q.extrapolate_first_order()
+        elif domain.order == 2:
+            Q.extrapolate_second_order()
+        else:
+            raise 'Unknown order'
+
+    # These quantities have been extrapolated.
+    # We may need to ensure top is never 0.
+    u_V  = velocity.vertex_values
+    w_V  = stage.vertex_values
+    h_V  = height.vertex_values
+    b_V  = width.vertex_values
+    t_V  = top.vertex_values
+        
+    # These need to be reconstructed
+    a_V  = area.vertex_values
+    d_V  = discharge.vertex_values
+    z_V  = bed.vertex_values
+
+    # State is true or false so should never be extrapolated.
+    s_V = state.vertex_values
+
+    # Reconstruct bed from stage and height.
+    z_V[:] = w_V-h_V
+
+    # Now reconstruct our conserved quantities from the above
+    # reconstructed quantities.
     a_V[:] = b_V*h_V
     d_V[:] = u_V*a_V
 
-    
     return
 

trunk/anuga_work/development/2010-projects/anuga_1d/pipe/pipe_domain_ext.c

@@ -21 +21 @@
 
 // Hard coded choices of computational methods
-#define flux_formula0 flux_formula_press0
-#define flux_formula1 flux_formula_press1
+//#define flux_formula0 flux_formula_press0
+//#define flux_formula1 flux_formula_press1
+#define flux_formula0 flux_formula_fs0
+#define flux_formula1 flux_formula_fs1
 #define compute_flux compute_flux_godunov
 #define forcingterm_formula_press forcingterm_formula_press_naive
 #define forcingterm_formula_fs forcingterm_formula_fs_naive
-#define compute_forcingterm compute_forcingterm_press
+//#define compute_forcingterm compute_forcingterm_press
+#define compute_forcingterm compute_forcingterm_fs
 
 // Some formulae
@@ -91 +94 @@
 
     edgeflux[1] = smax*flux_formula1(q_m, g) - smin*flux_formula1(q_p, g);
-    edgeflux[1] += smax*smin * (q_p[4]*a_p - q_m[4]*a_m);
+    //edgeflux[1] += smax*smin * (q_p[4]*a_p - q_m[4]*a_m);
+    edgeflux[1] += smax*smin * (q_p[1] - q_m[1]);
     edgeflux[1] /= denom;
     }

trunk/anuga_work/development/2010-projects/anuga_1d/pipe/test_pipe.py

@@ -12 +12 @@
 print "Pipe Variable Width Only Test"
 
+z_infty = 10.0       ## max equilibrium water depth at lowest point.
+L_x = 2500.0         ## width of channel
+A0 = 0.5*L_x                  ## determines amplitudes of oscillations
+omega = sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation
+
 # Define functions for initial quantities
 def initial_stage(x):
-    z_infty = 10.0       ## max equilibrium water depth at lowest point.
-    L_x = 2500.0         ## width of channel
-    A0 = 0.5*L_x                  ## determines amplitudes of oscillations
-    omega = sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation
     t=0.0
     y = numpy.zeros(len(x),numpy.float)
     for i in range(len(x)):
-        #y[i] = z_infty+2*A0*z_infty/L_x*cos(omega*t)*(x[i]/L_x-0.5*A0/(L_x)*cos(omega*t))
-        y[i] = 12.0
+        y[i] = z_infty+2*A0*z_infty/L_x*cos(omega*t)*(x[i]/L_x-0.5*A0/(L_x)*cos(omega*t))
+        #y[i] = 12.0
     return y
 
 def bed(x):
-    N = len(x)
-    z_infty = 10.0
-    z = numpy.zeros(N,numpy.float)
-    L_x = 2500.0
-    A0 = 0.5*L_x
-    omega = sqrt(2*g*z_infty)/L_x
-    for i in range(N):
+    z = numpy.zeros(len(x),numpy.float)
+    for i in range(len(x)):
         z[i] = z_infty*(x[i]**2/L_x**2)
     return z
+
+def width(x):
+    return 1.0
+
+def top(x):
+    return 4.0
 
 def initial_area(x):
@@ -40 +42 @@
     for i in range(len(x)):
         y[i]=initial_stage([x[i]])-bed([x[i]])
-
     return y
-
-def width(x):
-    return 1
-
-def top(x):
-    return 4
  
 def initial_state(x):
@@ -57 +52 @@
 
 # Set final time and yield time for simulation
-finaltime = 10.0
-yieldstep = 1.0
+finaltime = 500.0
+yieldstep = 10.0
 
 # Length of channel (m)