Changeset 8195

Timestamp:

Aug 10, 2011, 8:05:33 PM (14 years ago)

Author:

paul

Message:

channel_domain now has two distribute_to_vertices_and_edges functions available using different reconstruction methods

File:

• trunk/anuga_work/development/2010-projects/anuga_1d/channel/channel_domain.py (modified) (5 diffs)

trunk/anuga_work/development/2010-projects/anuga_1d/channel/channel_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)
         
 
@@ -164 +165 @@
 
 #-----------------------------------------------------------------------
-# Distribute to verticies with stage reconstructed and then extrapolated
+# Distribute to verticies with stage, velocity and channel 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
 
@@ -197 +192 @@
     w_C   = stage.centroid_values
 
-    #if domain.setstageflag:
-    #    a_C[:,] = (w_C-z_C)*b_C
-    #    domain.setstageflag = False
-
-    #if domain.discontinousb:
-     #   width.extrapolate_second_order()
-        
-
-    h0 = 1.0e-10
-    #Calculate height (and fix negatives)better be non-negative!
-    h_C[:]  = numpy.where(a_C > 0.0, a_C/b_C, 0.0)
-    u_C[:]  = numpy.where(a_C > 0.0, d_C/a_C, 0.0)
+    # Calculate height, velocity and stage.
+    # Here we assume the conserved quantities and the channel geometry
+    # (i.e. bed 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
-    
-
-    bed.extrapolate_second_order()
-    
-    for name in ['velocity', 'stage', 'width']:
+
+    # Extrapolate velocity and stage as well as channel geometry.
+    for name in ['velocity', 'stage', 'elevation', 'width']:
         Q = domain.quantities[name]
         if domain.order == 1:
@@ -234 +221 @@
     h_V  = height.vertex_values
     d_V  = discharge.vertex_values
-    
-
-
+
+    # Calculate height and fix up negatives. The main idea here is
+    # fix up the wet/dry interface.
     h_V[:,:] = w_V - z_V
 
@@ -253 +240 @@
 
 
+    # 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 edge values
-    w_V[:]  = z_V+h_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']
+    stage     = domain.quantities['stage']
+
+    #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
+    w_C   = stage.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']:
+        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.
+    u_V  = velocity.vertex_values
+    w_V  = stage.vertex_values
+    h_V  = height.vertex_values
+    b_V  = width.vertex_values
+
+    # These need to be reconstructed
+    a_V  = area.vertex_values
+    d_V  = discharge.vertex_values
+    z_V  = bed.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
 
-
-    print [ h_V, u_V ]
     return