Changeset 5727

Timestamp:

Sep 3, 2008, 10:28:15 PM (16 years ago)

Author:

steve

Message:

Setting up some more unit tests for shallow_water_domain

Location:

anuga_work/development/anuga_1d

Files:

• comp_flux_ext.c (modified) (2 diffs)
• dam_h_elevation.py (modified) (1 diff)
• shallow_water_domain.py (modified) (1 diff)
• shallow_water_domain_new.py (modified) (1 diff)
• test_shallow_water.py (modified) (9 diffs)

anuga_work/development/anuga_1d/comp_flux_ext.c

@@ -19 +19 @@
 /*      z=(a<b)?a:b; */
 /*      return z;} */
-
-<<<<<<< .mine
-
-
-
-
-=======
 
 
@@ -54 +47 @@
 
 
->>>>>>> .r5723
+
 //Innermost flux function (using w=z+h)
 int _flux_function(double *q_left, double *q_right,

anuga_work/development/anuga_1d/dam_h_elevation.py

@@ -1 +1 @@
 import os
 from math import sqrt, sin, cos, pi, exp
-from shallow_water_domain_new import *
+from shallow_water_domain import *
 from Numeric import zeros, Float
 #from analytic_dam_sudi import AnalyticDam

anuga_work/development/anuga_1d/shallow_water_domain.py

@@ -43 +43 @@
 """
 
-#from domain import *
-#from domain_order2 import *
+
 from domain import *
 Generic_Domain = Domain #Rename

anuga_work/development/anuga_1d/shallow_water_domain_new.py

@@ -154 +154 @@
         #Call correct module function
         #(either from this module or C-extension)
-        compute_fluxes_C_long(self) #compute_fluxes_C_wellbalanced(self)   #compute_fluxes_C_long(self)                 #compute_fluxes(self)    #compute_fluxes_C_short(self) 
+        compute_fluxes_C_short(self) #compute_fluxes_C_wellbalanced(self) #compute_fluxes_C_short(self)  #compute_fluxes_C_long(self)
 
     def compute_timestep(self):

anuga_work/development/anuga_1d/test_shallow_water.py

@@ -6 +6 @@
 
 from shallow_water_domain import *
+from shallow_water_domain import flux_function as domain_flux_function
+
 from Numeric import allclose, array, ones, Float
 
@@ -32 +34 @@
         domain.set_boundary({'exterior' : Reflective_boundary(domain)})
 
-        stage_ud, xmom_ud = compute_fluxes_python(domain)
+        stage_ud, xmom_ud = local_compute_fluxes(domain)
         
         domain.compute_fluxes()
 
-        #print doamin.quantities['xmomentum'].explicit_update
-        #print compute_fluxes_python(domain)
+        #print domain.quantities['stage'].explicit_update
+        #print domain.quantities['xmomentum'].explicit_update
+        print stage_ud
 
         assert allclose( domain.quantities['stage'].explicit_update, stage_ud )
@@ -51 +54 @@
 
         #This assumes h0 = 1.0e-3!!
-        edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr)
+        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
         #print maxspeed
         #print edgeflux
         
-        assert allclose(array([1.998002, 8.89201198],Float), edgeflux)
-        assert allclose(5.1284971665, maxspeed)
+        assert allclose(array([2.0, 8.9],Float), edgeflux, rtol=1.0e-005)
+        assert allclose(5.1305, maxspeed, rtol=1.0e-005)
 
         normal = -1.0
@@ -64 +67 @@
         zr = 0.0
 
-        edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr)
+        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
 
 
@@ -70 +73 @@
         #print edgeflux        
         
-        assert allclose(array([-1.998002, -8.89201198],Float), edgeflux)
-        assert allclose(5.1284971665, maxspeed)
+        assert allclose(array([-2.0, -8.9],Float), edgeflux, rtol=1.0e-005)
+        assert allclose(5.1305, maxspeed, rtol=1.0e-005)
 
     def test_domain_flux_function(self):
@@ -80 +83 @@
         zr = 0.0
 
-        edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr)
+        edgeflux, maxspeed = local_flux_function(normal, ql,qr,zl,zr)
 
         #print edgeflux
 
-        from shallow_water_domain import flux_function as domain_flux_function
+
 
         domainedgeflux, domainmaxspeed = domain_flux_function(normal, ql,qr,zl,zr)
@@ -93 +96 @@
         assert allclose(domainmaxspeed, maxspeed)
 
+    def test_gravity(self):
+        """
+        Compare shallow_water_domain gravity calculation
+        """
+                
+        def slope_one(x):
+            return x
+                        
+        domain = Domain(self.points)
+        domain.set_quantity('stage',4.0)
+        domain.set_quantity('elevation',slope_one)
+        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+        gravity(domain)
+        
+        #print domain.quantities['stage'].vertex_values 
+        #print domain.quantities['elevation'].vertex_values 
+        #print domain.quantities['xmomentum'].explicit_update       
+
+        assert allclose( array([-34.3, -24.5, -14.7], Float), domain.quantities['xmomentum'].explicit_update )
+
+
+    def test_evolve(self):
+        """
+        Compare shallow_water_domain gravity calculation
+        """
+                
+        def slope_one(x):
+            return x
+                        
+        domain = Domain(self.points)
+        domain.set_quantity('stage',4.0)
+        domain.set_quantity('elevation',slope_one)
+        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+        yieldstep=0.01
+        finaltime=0.01
+
+        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
+            domain.write_time()
+        
+        print domain.quantities['stage'].vertex_values 
+        print domain.quantities['elevation'].vertex_values 
+        print domain.quantities['xmomentum'].vertex_values
+
+
+        print domain.quantities['stage'].centroid_values 
+        print domain.quantities['elevation'].centroid_values 
+        print domain.quantities['xmomentum'].centroid_values    
+
+        assert allclose( array([-34.3, -24.5, -14.7], Float), domain.quantities['xmomentum'].explicit_update )
 
 
 #==============================================================================
 
-def compute_fluxes_python(domain):
+def local_compute_fluxes(domain):
     """Compute all fluxes and the timestep suitable for all volumes
     in domain.
@@ -122 +176 @@
     N = domain.number_of_elements
 
-    tmp0 = zeros((N,),Float)
-    tmp1 = zeros((N,),Float)
+    tmp0 = zeros(N,Float)
+    tmp1 = zeros(N,Float)
 
     #Shortcuts
@@ -222 +276 @@
 
 
-def flux_function(normal, ql, qr, zl, zr):
+def local_flux_function(normal, ql, qr, zl, zr):
     """Compute fluxes between volumes for the shallow water wave equation
     cast in terms of w = h+z using the 'central scheme' as described in