Changeset 7452 for anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

Timestamp:

Sep 1, 2009, 10:29:19 AM (15 years ago)

Author:

steve

Message:

Committing latest parallel code

File:

• anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py (modified) (1 diff)

anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -1644 +1644 @@
                                                              verbose=False)
                 assert num.allclose(Et, w + 0.5*u*u/g)
+
+
+    def test_cross_section_class(self):
+        """test_cross_section_class(self):
+
+        Test that the total and specific energy through a cross section can be
+        correctly obtained at run-time from the ANUGA cross section class.
+
+        This test creates a flat bed with a known flow through it, creates a cross
+        section and tests that the correct flow and energies are calculated
+
+        The specifics are
+        e = -1 m
+        u = 2 m/s
+        h = 2 m
+        w = 3 m (width of channel)
+
+        q = u*h*w = 12 m^3/s
+
+        This run tries it with georeferencing and with elevation = -1
+        """
+
+        import time, os
+        from Scientific.IO.NetCDF import NetCDFFile
+        from mesh_factory import rectangular
+
+        # Create basic mesh (20m x 3m)
+        width = 3
+        length = 20
+        t_end = 1
+        points, vertices, boundary = rectangular(length, width, length, width)
+
+        # Create shallow water domain
+        domain = Domain(points, vertices, boundary,
+                        geo_reference=Geo_reference(56, 308500, 6189000))
+
+        domain.default_order = 2
+        domain.set_quantities_to_be_stored(None)
+
+        e = -1.0
+        w = 1.0
+        h = w-e
+        u = 2.0
+        uh = u*h
+
+        Br = Reflective_boundary(domain)       # Side walls
+        Bd = Dirichlet_boundary([w, uh, 0])    # 2 m/s across the 3 m inlet:
+
+        # Initial conditions
+        domain.set_quantity('elevation', e)
+        domain.set_quantity('stage', w)
+        domain.set_quantity('xmomentum', uh)
+        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
+
+        # Interpolation points down the middle
+        I = [[0, width/2.],
+             [length/2., width/2.],
+             [length, width/2.]]
+        interpolation_points = domain.geo_reference.get_absolute(I)
+
+        # Shortcuts to quantites
+        stage = domain.get_quantity('stage')
+        xmomentum = domain.get_quantity('xmomentum')
+        ymomentum = domain.get_quantity('ymomentum')
+
+
+        # Create some cross sections
+        cross_sections = []
+        for i in range(5):
+            x = length/2. + i*0.23674563    # Arbitrary
+            polyline = [[x, 0], [x, width]]
+
+            polyline = domain.geo_reference.get_absolute(polyline)
+        
+            cross_sections.append(Cross_section(domain,polyline))
+
+
+ 
+        for t in domain.evolve(yieldstep=0.1, finaltime=t_end):
+            # Check that quantities are they should be in the interior
+            w_t = stage.get_values(interpolation_points)
+            uh_t = xmomentum.get_values(interpolation_points)
+            vh_t = ymomentum.get_values(interpolation_points)
+
+            assert num.allclose(w_t, w)
+            assert num.allclose(uh_t, uh)
+            assert num.allclose(vh_t, 0.0, atol=1.0e-6)
+
+
+            # Check flows and energies through the middle
+            for cross_section in cross_sections:
+                
+                Q = cross_section.get_flow_through_cross_section()
+
+                assert num.allclose(Q, uh*width)
+
+                Es = cross_section.get_energy_through_cross_section(kind='specific')
+
+                assert num.allclose(Es, h + 0.5*u*u/g)
+
+                Et = cross_section.get_energy_through_cross_section(kind='total')
+                
+                assert num.allclose(Et, w + 0.5*u*u/g)
+
+
+
+
+
 
     def test_another_runup_example(self):