Python Domain.set_minimum_storable_height方法代码示例

# 需要导入模块: from anuga.shallow_water.shallow_water_domain import Domain [as 别名]
# 或者: from anuga.shallow_water.shallow_water_domain.Domain import set_minimum_storable_height [as 别名]
    def test_get_flow_through_cross_section_stored_uniquely(self):
        """test_get_flow_through_cross_section_stored_uniquely(self):

        Test that the total flow through a cross section can be
        correctly obtained from an sww file.
        
        This test creates a flat bed with a known flow through it and tests
        that the function correctly returns the expected flow.

        The specifics are
        u = 2 m/s
        h = 1 m
        w = 3 m (width of channel)

        q = u*h*w = 6 m^3/s
       
        
        """

        import time, os
        from anuga.file.netcdf import NetCDFFile

        # Setup
        #from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular

        # Create basic mesh (20m x 3m)
        width = 3
        length = 20
        t_end = 3
        points, vertices, boundary = rectangular(length, width,
                                                 length, width)

        # Create shallow water domain
        domain = Domain(points, vertices, boundary)
        domain.default_order = 2
        domain.set_minimum_storable_height(0.01)

        domain.set_name('flowtest_uniquely')
        swwfile = domain.get_name() + '.sww'

        domain.set_store_vertices_uniquely()
        
        domain.set_datadir('.')
        domain.format = 'sww'
        domain.smooth = True

        h = 1.0
        u = 2.0
        uh = u*h

        Br = Reflective_boundary(domain)     # Side walls
        Bd = Dirichlet_boundary([h, uh, 0])  # 2 m/s across the 3 m inlet: 


        
        domain.set_quantity('elevation', 0.0)
        domain.set_quantity('stage', h)
        domain.set_quantity('xmomentum', uh)
        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})

        for t in domain.evolve(yieldstep=1, finaltime = t_end):
            pass

        # Check that momentum is as it should be in the interior

        I = [[0, width/2.],
             [length/2., width/2.],
             [length, width/2.]]
        
        f = file_function(swwfile,
                          quantities=['stage', 'xmomentum', 'ymomentum'],
                          interpolation_points=I,
                          verbose=False)
        for t in range(t_end+1):
            for i in range(3):
                assert num.allclose(f(t, i), [1, 2, 0], atol=1.0e-6)
            

        # Check flows through the middle
        for i in range(5):
            x = length/2. + i*0.23674563 # Arbitrary
            cross_section = [[x, 0], [x, width]]
            time, Q = get_flow_through_cross_section(swwfile,
                                                     cross_section,
                                                     verbose=False)

            assert num.allclose(Q, uh*width)


       
        # Try the same with partial lines
        x = length/2.
        for i in range(5):
            start_point = [length/2., i*width/5.]
            #print start_point
                            
            cross_section = [start_point, [length/2., width]]
            time, Q = get_flow_through_cross_section(swwfile,
                                                     cross_section,
                                                     verbose=False)
#.........这里部分代码省略.........

# 需要导入模块: from anuga.shallow_water.shallow_water_domain import Domain [as 别名]
# 或者: from anuga.shallow_water.shallow_water_domain.Domain import set_minimum_storable_height [as 别名]
    def test_get_energy_through_cross_section(self):
        """test_get_energy_through_cross_section(self):

        Test that the specific and total energy through a cross section can be
        correctly obtained from an sww file.
        
        This test creates a flat bed with a known flow through it and tests
        that the function correctly returns the expected energies.

        The specifics are
        u = 2 m/s
        h = 1 m
        w = 3 m (width of channel)

        q = u*h*w = 6 m^3/s
        Es = h + 0.5*v*v/g  # Specific energy head [m]
        Et = w + 0.5*v*v/g  # Total energy head [m]        


        This test uses georeferencing
        
        """

        import time, os
        from anuga.file.netcdf import NetCDFFile

        # Setup
        #from anuga.abstract_2d_finite_volumes.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_minimum_storable_height(0.01)

        domain.set_name('flowtest')
        swwfile = domain.get_name() + '.sww'

        domain.set_datadir('.')
        domain.format = 'sww'
        domain.smooth = True

        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: 

        
        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})

        for t in domain.evolve(yieldstep=1, finaltime = t_end):
            pass

        # Check that momentum is as it should be in the interior

        I = [[0, width/2.],
             [length/2., width/2.],
             [length, width/2.]]
        
        I = domain.geo_reference.get_absolute(I)
        f = file_function(swwfile,
                          quantities=['stage', 'xmomentum', 'ymomentum'],
                          interpolation_points=I,
                          verbose=False)

        for t in range(t_end+1):
            for i in range(3):
                #print i, t, f(t, i)
                assert num.allclose(f(t, i), [w, uh, 0], atol=1.0e-6)
            

        # Check energies through the middle
        for i in range(5):
            x = length/2. + i*0.23674563 # Arbitrary
            cross_section = [[x, 0], [x, width]]

            cross_section = domain.geo_reference.get_absolute(cross_section)            
            
            time, Es = get_energy_through_cross_section(swwfile,
                                                       cross_section,
                                                       kind='specific',
                                                       verbose=False)
            assert num.allclose(Es, h + 0.5*u*u/g)
            
            time, Et = get_energy_through_cross_section(swwfile,
#.........这里部分代码省略.........

# 需要导入模块: from anuga.shallow_water.shallow_water_domain import Domain [as 别名]
# 或者: from anuga.shallow_water.shallow_water_domain.Domain import set_minimum_storable_height [as 别名]
    def test_get_maximum_inundation_de0(self):
        """Test that sww information can be converted correctly to maximum
        runup elevation and location (without and with georeferencing)

        This test creates a slope and a runup which is maximal (~11m) at around 10s
        and levels out to the boundary condition (1m) at about 30s.
        """

        import time, os
        from anuga.file.netcdf import NetCDFFile

        #Setup

        #from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular

        # Create basic mesh (100m x 100m)
        points, vertices, boundary = rectangular(20, 5, 100, 50)

        # Create shallow water domain
        domain = Domain(points, vertices, boundary)
        domain.default_order = 2
        domain.set_minimum_storable_height(0.01)

        filename = 'runup_test_3'
        domain.set_name(filename)
        swwfile = domain.get_name() + '.sww'

        domain.set_datadir('.')
        domain.format = 'sww'
        domain.smooth = True

        # FIXME (Ole): Backwards compatibility
        # Look at sww file and see what happens when
        # domain.tight_slope_limiters = 1
        domain.tight_slope_limiters = 0
        domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)        
        
        Br = Reflective_boundary(domain)
        Bd = Dirichlet_boundary([1.0,0,0])


        #---------- First run without geo referencing
        
        domain.set_quantity('elevation', lambda x,y: -0.2*x + 14) # Slope
        domain.set_quantity('stage', -6)
        domain.set_boundary( {'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})

        for t in domain.evolve(yieldstep=1, finaltime = 50):
            pass


        # Check maximal runup
        runup = get_maximum_inundation_elevation(swwfile)
        location = get_maximum_inundation_location(swwfile)
        #print 'Runup, location', runup, location
        assert num.allclose(runup, 4.66666666667)
        assert num.allclose(location[0], 46.666668) 
               
        # Check final runup
        runup = get_maximum_inundation_elevation(swwfile, time_interval=[45,50])
        location = get_maximum_inundation_location(swwfile, time_interval=[45,50])
        #print 'Runup, location:',runup, location

        assert num.allclose(runup, 3.81481488546)
        assert num.allclose(location[0], 51.666668)

        # Check runup restricted to a polygon
        p = [[50,1], [99,1], [99,49], [50,49]]
        runup = get_maximum_inundation_elevation(swwfile, polygon=p)
        location = get_maximum_inundation_location(swwfile, polygon=p)
        #print runup, location

        assert num.allclose(runup, 3.81481488546) 
        assert num.allclose(location[0], 51.6666666)                

        # Check that mimimum_storable_height works
        fid = NetCDFFile(swwfile, netcdf_mode_r) # Open existing file
        
        stage = fid.variables['stage_c'][:]
        z = fid.variables['elevation_c'][:]
        xmomentum = fid.variables['xmomentum_c'][:]
        ymomentum = fid.variables['ymomentum_c'][:]
        
        for i in range(stage.shape[0]):
            h = stage[i]-z # depth vector at time step i
            
            # Check every node location
            for j in range(stage.shape[1]):
                # Depth being either exactly zero implies
                # momentum being zero.
                # Or else depth must be greater than or equal to
                # the minimal storable height
                if h[j] == 0.0:
                    assert xmomentum[i,j] == 0.0
                    assert ymomentum[i,j] == 0.0                
                else:
                    assert h[j] >= 0.0
        
        fid.close()

#.........这里部分代码省略.........