Changeset 8382

Timestamp:

Apr 3, 2012, 8:26:00 PM (13 years ago)

Author:

steve

Message:

can run anuga with various flux computations (2 well balanced schemes)

Location:

trunk

Files:

• anuga_core/source/anuga/shallow_water/shallow_water_ext.c (modified) (4 diffs)
• anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py (modified) (3 diffs)
• anuga_work/shallow_water_balanced_steve/run_sw_merimbula.py (modified) (2 diffs)

trunk/anuga_core/source/anuga/shallow_water/shallow_water_ext.c

@@ -2913 +2913 @@
 
     // Local variables
-    double max_speed, length, inv_area, zl, zr;
+    double max_speed, length, inv_area;
     double timestep = 1.0e30;
     double h0 = D->H0*D->H0; // This ensures a good balance when h approaches H0.
@@ -2927 +2927 @@
     double hl, hl1, hl2;
     double hr, hr1, hr2;
+    double zl, zl1, zl2;
+    double zr, zr1, zr2;
+    double wr;
 
     // Workspace (making them static actually made function slightly slower (Ole))
@@ -2966 +2969 @@
 
             zl = D->bed_edge_values[k3i];
+
+            zl1 = D->bed_vertex_values[k3i1];
+            zl2 = D->bed_vertex_values[k3i2];
+
+
             hl = D->stage_edge_values[k3i] - zl;
 
-            hl1 = D->stage_vertex_values[k3i1] - D->bed_vertex_values[k3i1];
-            hl2 = D->stage_vertex_values[k3i2] - D->bed_vertex_values[k3i2];
+            hl1 = D->stage_vertex_values[k3i1] - zl1;
+            hl2 = D->stage_vertex_values[k3i2] - zl2;
+
 
             // Get right hand side values either from neighbouring triangle
@@ -2982 +2991 @@
                 qr[2] = D->ymom_boundary_values[m];
 
-                zr = zl; // Extend bed elevation to boundary and assume flat
-                hr = D->stage_boundary_values[m] - zr;
-                hr1 = hr;
-                hr2 = hr;
+                zr = zl; // Extend bed elevation to boundary and assume flat stage
+
+                wr = D->stage_boundary_values[m];
+
+                hr = wr - zr;
+                hr1 = wr-zl2;
+                hr2 = wr-zl1;
+
 
             }

trunk/anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -2227 +2227 @@
         assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0)
 
-    def test_gravity_3(self):
-        #Showing standard gravity term is not well balanced
+
+
+    def test_well_balanced_flux_1(self):
+        #testing that compute_fluxes_method = wb_1 is well balance
 
         from anuga.config import g
+        from anuga import Reflective_boundary
 
         a = [0.0, 0.0]
@@ -2256 +2259 @@
         domain.set_quantity('stage', stage)
 
+
+
         for name in domain.conserved_quantities:
             assert num.allclose(domain.quantities[name].explicit_update, 0)
             assert num.allclose(domain.quantities[name].semi_implicit_update, 0)
 
-        domain.set_compute_fluxes_method('original')
-        domain.compute_forcing_terms()
+
+        Br = Reflective_boundary(domain)      # Solid reflective wall
+        domain.set_boundary({'exterior' :Br})
+        domain.update_boundary()
+
+        domain.set_compute_fluxes_method('wb_1')
+        domain.compute_fluxes()
 
 
         print domain.quantities['xmomentum'].explicit_update
-        #print domain.quantities['ymomentum'].explicit_update
+        print domain.quantities['stage'].vertex_values
+        print domain.quantities['elevation'].vertex_values
+        print domain.quantities['ymomentum'].explicit_update
 
 
         assert num.allclose(domain.quantities['stage'].explicit_update, 0)
-        assert num.allclose(domain.quantities['xmomentum'].explicit_update,
-                            [-382.2, -323.4, -205.8, -382.2])
-        assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0)
-
+        assert num.allclose(domain.quantities['xmomentum'].explicit_update, 0.0)
+        assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0.0)
+
+
+    def test_well_balanced_flux_2(self):
+        #esting that compute_fluxes_method = wb_2 is well balance
+
+        from anuga.config import g
+        from anuga import Reflective_boundary
+
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #             bac,     bce,     ecf,     dbe
+        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        domain = Domain(points, vertices)
+
+        #Set up for a gradient of (3,0) at mid triangle (bce)
+        def slope(x, y):
+            return 3*x
+
+        h = 15
+        def stage(x, y):
+            return h
+
+        domain.set_quantity('elevation', slope)
+        domain.set_quantity('stage', stage)
+
+
+
+        for name in domain.conserved_quantities:
+            assert num.allclose(domain.quantities[name].explicit_update, 0)
+            assert num.allclose(domain.quantities[name].semi_implicit_update, 0)
+
+
+        Br = Reflective_boundary(domain)      # Solid reflective wall
+        domain.set_boundary({'exterior' :Br})
+        domain.update_boundary()
+
+        domain.set_compute_fluxes_method('wb_2')
+        domain.compute_fluxes()
+
+
+        print domain.quantities['xmomentum'].explicit_update
+        print domain.quantities['stage'].vertex_values
+        print domain.quantities['elevation'].vertex_values
+        print domain.quantities['ymomentum'].explicit_update
+
+
+        assert num.allclose(domain.quantities['stage'].explicit_update, 0)
+        assert num.allclose(domain.quantities['xmomentum'].explicit_update, 0.0)
+        assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0.0)
 
     def test_gravity_wb(self):
@@ -7585 +7651 @@
 if __name__ == "__main__":
     #suite = unittest.makeSuite(Test_Shallow_Water, 'test_rainfall_forcing_with_evolve')
-    suite = unittest.makeSuite(Test_Shallow_Water, 'test')
+    suite = unittest.makeSuite(Test_Shallow_Water, 'test_well')
     runner = unittest.TextTestRunner(verbosity=1)
     runner.run(suite)

trunk/anuga_work/shallow_water_balanced_steve/run_sw_merimbula.py

@@ -49 +49 @@
 #mesh_filename = "test-100.tsh" ; x0 = 0.25 ; x1 = 0.5
 #mesh_filename = "test-20.tsh" ; x0 = 250.0 ; x1 = 350.0
-yieldstep = 1 #20
-finaltime = 10 #200
+yieldstep = 20
+finaltime = 200
 verbose = True
 
@@ -97 +97 @@
 #domain.set_CFL(0.7)
 #domain.set_beta(1.5)
-domain.set_compute_fluxes_method('wb_2')
+domain.set_compute_fluxes_method('wb_1')
 domain.set_name('merimbula')