Changeset 229 for inundation/ga/storm_surge/pyvolution/shallow_water.py

Timestamp:

Aug 26, 2004, 10:14:19 PM (20 years ago)

Author:

ole

Message:

Finally quashed bug where ymomentum wasn't updated in
second order limiter

File:

• inundation/ga/storm_surge/pyvolution/shallow_water.py (modified) (6 diffs)

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -273 +273 @@
 def distribute_to_vertices_and_edges(domain):
 
-    #print 'Calling distrib within SW'
+    ##print 'Calling distrib within SW'
     if domain.order == 1:
+        #FIXME: This should be cleaned up, but we try to follow 
+        #pyvolution 2 strictly for now
+        protect_against_negative_heights_centroid(domain)        
         protect_against_infinitesimal_heights_centroid(domain)        
         extrapolate_first_order(domain)
     elif domain.order == 2:
-        protect_against_negative_heights_centroid(domain)        
+        #protect_against_negative_heights_centroid(domain)        
         extrapolate_second_order(domain)        
     else:
@@ -323 +326 @@
     """
 
-    #FIXME: USed only in second order
+    
     #Water levels at centroids
     wc = domain.quantities['level'].centroid_values
@@ -428 +431 @@
         #Update water levels
         for i in range(3):
-            #wv[k,i] = zv[k,i] + hc[k] + alpha*(zc[k]-zv[k,i])
-            wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k])
-
-
+            #FIXME: Use the original first-order one first, then switch
+            wv[k,i] = zv[k,i] + hc[k] + alpha*(zc[k]-zv[k,i])
+            #wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k])
+
+        #FIXME: What about alpha weighting of momentum??
+            
         
 def extrapolate_second_order(domain):
@@ -462 +467 @@
         Q = domain.quantities[name]
         Q.extrapolate_second_order()
+        
+    #print 'y1', Q.vertex_values[1,:]   #OK
+    
+    #FIXME - like pyvolution 2 .......................
+    protect_against_negative_heights_centroid(domain)   
+
+    #print 'y1', Q.vertex_values[1,:]   #OK    
+    
+    for name in domain.conserved_quantities:
+        Q = domain.quantities[name]     
         Q.limit()
 
   
+    #print 'y1', Q.vertex_values[1,:]   #OK    
+      
     #Water levels at centroids
     wc = domain.quantities['level'].centroid_values
@@ -507 +524 @@
         #first order scheme and alpha==1 means using the limited
         #second order scheme for the stage w
-        #If hmin < 0 then alpha=0 reverrting to first order.
+        #If hmin < 0 then alpha=0 reverting to first order.
       
         if z_range > 0.0:
@@ -513 +530 @@
         else: 
             alpha = 1.0
-
+   
+        ##if k==1: print 'alpha', alpha #OK
+        
         #Update water levels
         #  (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) = 
-        #   zvi + hc + alpha*(hvi - hc)         
-        for i in range(3):
-            wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k])
-            
-
-        #Momentums at centroids
-        xmomc = domain.quantities['xmomentum'].centroid_values
-        ymomc = domain.quantities['ymomentum'].centroid_values        
-
-        #Momentums at vertices
-        xmomv = domain.quantities['xmomentum'].vertex_values
-        ymomv = domain.quantities['ymomentum'].vertex_values                
-
-        # Update momentum as a linear combination of 
-        # xmomc and ymomc (shallow) and momentum
-        # from extrapolator xmomv and ymomv (deep).
-        
-        xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:];            
-            
-
+        #   zvi + hc + alpha*(hvi - hc)
+        if alpha < 1:         
+            for i in range(3):
+                wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k])
+            
+
+            #Momentums at centroids
+            xmomc = domain.quantities['xmomentum'].centroid_values
+            ymomc = domain.quantities['ymomentum'].centroid_values        
+
+            #Momentums at vertices
+            xmomv = domain.quantities['xmomentum'].vertex_values
+            ymomv = domain.quantities['ymomentum'].vertex_values         
+
+            # Update momentum as a linear combination of 
+            # xmomc and ymomc (shallow) and momentum
+            # from extrapolator xmomv and ymomv (deep).
+        
+            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:];            
+            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:];                     
+            
+    #print 'y1', Q.vertex_values[1,:]   #OK    
+            
     #Finally, protect against infinitesimal heights and high speeds
     #Water depths at vertices