Changeset 5731 for anuga_work/development/anuga_1d/test_shallow_water.py

Timestamp:

Sep 4, 2008, 8:28:49 PM (15 years ago)

Author:

steve

Message:

Unit test for second order timestepping

File:

• anuga_work/development/anuga_1d/test_shallow_water.py (modified) (3 diffs)

anuga_work/development/anuga_1d/test_shallow_water.py

@@ -5 +5 @@
 
 
-from shallow_water_domain_new import *
-from shallow_water_domain_new import flux_function as domain_flux_function
-
-from Numeric import allclose, array, ones, Float, maximum
+from shallow_water_domain import *
+from shallow_water_domain import flux_function as domain_flux_function
+
+from Numeric import allclose, array, ones, Float, maximum, zeros
 
 
@@ -119 +119 @@
 
 
-    def test_evolve(self):
-        """
-        Compare shallow_water_domain gravity calculation
+    def test_evolve_first_order(self):
+        """
+        Compare still lake solution for various versions of shallow_water_domain
         """
                 
@@ -132 +132 @@
         domain.set_boundary({'exterior' : Reflective_boundary(domain)})
 
+        domain.default_order = 1
+        domain.default_time_order = 1
         yieldstep=10.0
         finaltime=10.0
 
         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 )
+            pass
+        
+##      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( 10.0*ones(10), domain.quantities['stage'].centroid_values )
+        assert allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
+
+
+    def test_evolve_euler_second_order_space(self):
+        """
+        Compare still lake solution for various versions of shallow_water_domain
+        """
+                
+        def slope_square(x):
+            return maximum(4.0-(x-5.0)*(x-5.0), 0.0)
+                        
+        domain = Domain(self.points2)
+        domain.set_quantity('stage',10.0)
+        domain.set_quantity('elevation',slope_square)
+        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+        domain.default_order = 2
+        domain.default_time_order = 1
+        yieldstep=10.0
+        finaltime=10.0
+
+        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
+            pass
+        
+        assert allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
+        assert allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
+
+    def test_evolve_euler_second_order_space_time(self):
+        """
+        Compare still lake solution for various versions of shallow_water_domain
+        """
+                
+        def slope_square(x):
+            return maximum(4.0-(x-5.0)*(x-5.0), 0.0)
+                        
+        domain = Domain(self.points2)
+        domain.set_quantity('stage',10.0)
+        domain.set_quantity('elevation',slope_square)
+        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
+
+        domain.default_order = 2
+        domain.default_time_order = 2
+        yieldstep=10.0
+        finaltime=10.0
+
+        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
+            pass
+        
+        assert allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
+        assert allclose( zeros(10), domain.quantities['xmomentum'].centroid_values )
+