Changeset 7868 for anuga_work/development/2010-projects/anuga_1d/sww/sww_domain.py

Timestamp:

Jun 22, 2010, 5:30:32 PM (13 years ago)

Author:

steve

Message:

Added in some more c based limiters

File:

• anuga_work/development/2010-projects/anuga_1d/sww/sww_domain.py (modified) (5 diffs)

anuga_work/development/2010-projects/anuga_1d/sww/sww_domain.py

@@ -90 +90 @@
         self.quantities_to_be_stored = ['stage','xmomentum']
 
-        self.__doc__ = 'shallow_water_domain'
+        self.__doc__ = 'sww_domain'
 
 
@@ -178 +178 @@
             yield(t)
 
-    def initialise_storage(self):
-        """Create and initialise self.writer object for storing data.
-        Also, save x and bed elevation
-        """
-
-        import data_manager
-
-        #Initialise writer
-        self.writer = data_manager.get_dataobject(self, mode = 'w')
-
-        #Store vertices and connectivity
-        self.writer.store_connectivity()
-
-
-    def store_timestep(self, name):
-        """Store named quantity and time.
-
-        Precondition:
-           self.write has been initialised
-        """
-        self.writer.store_timestep(name)
+
+
 
 
@@ -253 +234 @@
     u_C   = Velocity.centroid_values
         
-    #print id(h_C)
-    #FIXME replace with numpy.where
-    for i in range(N):
-        h_C[i] = w_C[i] - z_C[i]                                                
-        if h_C[i] <= 0.0:
-            #print 'h_C[%d]= %15.5e\n' % (i,h_C[i])
-            h_C[i] = 1.0e-15
-            w_C[i] = z_C[i]
-            uh_C[i] = 0.0
-                
-    
-##    for i in range(len(h_C)):
-##      if h_C[i] < epsilon:
-##          u_C[i] = 0.0  #Could have been negative
-##          h_C[i] = 0.0
-##      else:
-
-    u_C[:,]  = uh_C/(h_C + h0/h_C)
+    #Calculate height (and fix negatives)better be non-negative!
+    h_C[:] = w_C - z_C
+    u_C[:]  = uh_C/(h_C + h0/h_C)
         
     for name in [ 'velocity', 'stage' ]:
@@ -284 +250 @@
             raise 'Unknown order'
 
+
     stage_V = domain.quantities['stage'].vertex_values                  
     bed_V   = domain.quantities['elevation'].vertex_values      
@@ -291 +258 @@
 
     h_V[:,:]    = stage_V - bed_V
-    for i in range(len(h_C)):
-        for j in range(2):
-            if h_V[i,j] < 0.0 :
-                #print 'h_V[%d,%d] = %f \n' % (i,j,h_V[i,j])                 
-                dh = h_V[i,j]
-                h_V[i,j] = 0.0
-                stage_V[i,j] = bed_V[i,j]
-                h_V[i,(j+1)%2] = h_V[i,(j+1)%2] + dh
-                stage_V[i,(j+1)%2] = stage_V[i,(j+1)%2] + dh
-                
+
+#    # protect from edge values going negative
+#    h_V[:,1] = numpy.where(h_V[:,0] < 0.0 , h_V[:,1]-h_V[:,0], h_V[:,1])
+#    h_V[:,0] = numpy.where(h_V[:,0] < 0.0 , 0.0, h_V[:,0])
+#
+#    h_V[:,0] = numpy.where(h_V[:,1] < 0.0 , h_V[:,0]-h_V[:,1], h_V[:,0])
+#    h_V[:,1] = numpy.where(h_V[:,1] < 0.0 , 0.0, h_V[:,1])
+#
+#
+#    stage_V[:,:] = bed_V + h_V
     xmom_V[:,:] = u_V * h_V