Changeset 5536

Timestamp:

Jul 18, 2008, 4:37:03 PM (15 years ago)

Author:

steve

Message:

Added unit test files for quantity and shallow water

Location:

anuga_work/development/anuga_1d

Files:

• domain.py (modified) (2 diffs)
• dry_dam_sudi.py (modified) (1 diff)
• quantity.py (modified) (7 diffs)
• shallow_water_domain.py (modified) (7 diffs)
• test_quantity.py (added)
• test_shallow_water.py (added)

anuga_work/development/anuga_1d/domain.py

@@ -34 +34 @@
         self.beta = 1.0
         self.limiter = "minmod_kurganov"
-        self.split = False
         self.wet_nodes = zeros((N,2), Int) # should this be here
 
@@ -125 +124 @@
         #FIXME: remove later - maybe OK, though....
         for name in self.conserved_quantities:
-            self.quantities[name] = Conserved_quantity(self)
+            self.quantities[name] = Quantity(self)
         for name in self.other_quantities:
             self.quantities[name] = Quantity(self)

anuga_work/development/anuga_1d/dry_dam_sudi.py

@@ -1 +1 @@
 import os
 from math import sqrt, pi
-from shallow_water_1d import *
+from shallow_water_domain import *
 from Numeric import allclose, array, zeros, ones, Float, take, sqrt
 from config import g, epsilon

anuga_work/development/anuga_1d/quantity.py

@@ -22 +22 @@
         #from domain import Domain
         #from domain_order2 import Domain
-        from domain_t2 import Domain
+        from domain import Domain
         from Numeric import array, zeros, Float
 
@@ -56 +56 @@
         #Intialise centroid values
         self.interpolate()
+
+
+        from Numeric import zeros, Float
+        
+        #Allocate space for boundary values
+        #L = len(domain.boundary)
+        self.boundary_values = zeros(2, Float) #assumes no parrellism
+
+        #Allocate space for updates of conserved quantities by
+        #flux calculations and forcing functions
+        
+        N = domain.number_of_elements
+        self.explicit_update = zeros(N, Float )
+        self.semi_implicit_update = zeros(N, Float )
+
+        self.gradients = zeros(N, Float)
+        self.qmax = zeros(self.centroid_values.shape, Float)
+        self.qmin = zeros(self.centroid_values.shape, Float)
+
+        self.beta = domain.beta        
 
     #Methods for operator overloading
@@ -381 +401 @@
         return integral
 
-class Conserved_quantity(Quantity):
-    """Class conserved quantity adds to Quantity:
-
-    storage and method for updating, and
-    methods for extrapolation from centropid to vertices inluding
-    gradients and limiters
-    """
-
-    def __init__(self, domain, vertex_values=None):
-        Quantity.__init__(self, domain, vertex_values)
-
-        from Numeric import zeros, Float
-        
-        #Allocate space for boundary values
-        #L = len(domain.boundary)
-        self.boundary_values = zeros(2, Float) #assumes no parrellism
-
-        #Allocate space for updates of conserved quantities by
-        #flux calculations and forcing functions
-        
-        N = domain.number_of_elements
-        self.explicit_update = zeros(N, Float )
-        self.semi_implicit_update = zeros(N, Float )
-
-        self.gradients = zeros(N, Float)
-        self.qmax = zeros(self.centroid_values.shape, Float)
-        self.qmin = zeros(self.centroid_values.shape, Float)
-
-        self.beta = domain.beta
 
     def update(self, timestep):
@@ -863 +854 @@
                 else:
                     qv[k,i] = qc[k]
+    
+
+class Conserved_quantity(Quantity):
+    """Class conserved quantity adds to Quantity:
+
+    storage and method for updating, and
+    methods for extrapolation from centropid to vertices inluding
+    gradients and limiters
+    """
+
+    def __init__(self, domain, vertex_values=None):
+        Quantity.__init__(self, domain, vertex_values)
+
+        print "Use Quantity instead of Conserved_quantity"
+
                     
 
@@ -884 +890 @@
 if __name__ == "__main__":
     #from domain import Domain
-    from shallow_water_1d import Domain     
+    from shallow_water_domain import Domain     
     from Numeric import arange
     
@@ -892 +898 @@
     D1 = Domain(points1)
 
-    Q1 = Conserved_quantity(D1, vertex_values)
+    Q1 = Quantity(D1, vertex_values)
 
     print Q1.vertex_values
@@ -941 +947 @@
     D2 = Domain(points2)
 
-    Q2 = Conserved_quantity(D2)
+    Q2 = Quantity(D2)
     Q2.set_values(fun,'vertices')
     Xc = Q2.domain.vertices

anuga_work/development/anuga_1d/shallow_water_domain.py

@@ -45 +45 @@
 #from domain import *
 #from domain_order2 import *
-from domain_t2 import *
+from domain import *
 Generic_Domain = Domain #Rename
 
@@ -388 +388 @@
 
 
-def flux_function_split(normal, ql, qr, zl, zr):
-    from config import g, epsilon
-    from math import sqrt
-    from Numeric import array
-
-    #print 'ql',ql
-
-    #Align momentums with x-axis
-    #q_left  = rotate(ql, normal, direction = 1)
-    #q_right = rotate(qr, normal, direction = 1)
-    q_left = ql
-    q_left[1] = q_left[1]*normal
-    q_right = qr
-    q_right[1] = q_right[1]*normal
-
-    #z = (zl+zr)/2 #Take average of field values
-    z = 0.5*(zl+zr) #Take average of field values
-
-    w_left  = q_left[0]   #w=h+z
-    h_left  = w_left-z
-    uh_left = q_left[1]
-
-    if h_left < epsilon:
-        u_left = 0.0  #Could have been negative
-        h_left = 0.0
-    else:
-        u_left  = uh_left/h_left
-
-
-    w_right  = q_right[0]  #w=h+z
-    h_right  = w_right-z
-    uh_right = q_right[1]
-
-
-    if h_right < epsilon:
-        u_right = 0.0  #Could have been negative
-        h_right = 0.0
-    else:
-        u_right  = uh_right/h_right
-
-    #vh_left  = q_left[2]
-    #vh_right = q_right[2]
-
-    #soundspeed_left  = sqrt(g*h_left)
-    #soundspeed_right = sqrt(g*h_right)
-
-    #Maximal wave speed
-    #s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
-    s_max = max(u_left, u_right, 0)
-    
-    #Minimal wave speed
-    #s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
-    s_min = min(u_left, u_right, 0)
-    
-    #Flux computation
-
-    #flux_left  = array([u_left*h_left,
-    #                    u_left*uh_left + 0.5*g*h_left*h_left])
-    #flux_right = array([u_right*h_right,
-    #                    u_right*uh_right + 0.5*g*h_right*h_right])
-    flux_left  = array([u_left*h_left,
-                        u_left*uh_left])# + 0.5*g*h_left*h_left])
-    flux_right = array([u_right*h_right,
-                        u_right*uh_right])# + 0.5*g*h_right*h_right])
-    
-    denom = s_max-s_min
-    if denom == 0.0:
-        edgeflux = array([0.0, 0.0])
-        max_speed = 0.0
-    else:
-        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
-        edgeflux += s_max*s_min*(q_right-q_left)/denom
-        
-        edgeflux[1] = edgeflux[1]*normal
-
-        max_speed = max(abs(s_max), abs(s_min))
-
-    return edgeflux, max_speed    
+  
 
 def compute_timestep(domain):
@@ -518 +441 @@
             normal = domain.normals[k, i]
 
-            if domain.split == False:
-                edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
-            elif domain.split == True:
-                edgeflux, max_speed = flux_function_split(normal, ql, qr, zl, zr)
+            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
+
             #Update optimal_timestep
             try:
@@ -612 +533 @@
                 #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
                 qr[0] = stage[n, m]
-                qr[1] =  xmom[n, m]
+                qr[1] = xmom[n, m]
                 zr = bed[n, m]
 
@@ -625 +546 @@
             
 
-            if domain.split == False:
-                edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
-            elif domain.split == True:
-                edgeflux, max_speed = flux_function_split(normal, ql, qr, zl, zr)
+            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
+
             #print 'edgeflux', edgeflux
 
@@ -646 +565 @@
 
         Stage.explicit_update[k] = flux[0]
-        Xmom.explicit_update[k] = flux[1]
+        Xmom.explicit_update[k] =  flux[1]
         #Ymom.explicit_update[k] = flux[2]
         #print "flux cell",k,flux[0]
@@ -1054 +973 @@
         
         #Update momentum (explicit update is reset to source values)
-        if domain.split == False:
-            xmom[k] += -g*bx*avg_h
-            #xmom[k] = -g*bx*avg_h
-            #stage[k] = 0.0
-        elif domain.split == True:
-            xmom[k] += -g*wx*avg_h
-            #xmom[k] = -g*wx*avg_h
-        #ymom[k] += -g*zy*avg_h
+        xmom[k] += -g*bx*avg_h
+        #xmom[k] = -g*bx*avg_h
+        #stage[k] = 0.0
+ 
  
 def manning_friction(domain):