Python Domain.evolve方法代码示例

# 需要导入模块: from anuga.shallow_water.shallow_water_domain import Domain [as 别名]
# 或者: from anuga.shallow_water.shallow_water_domain.Domain import evolve [as 别名]
    def test_merge_swwfiles(self):
        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular, \
                                                                    rectangular_cross
        from anuga.shallow_water.shallow_water_domain import Domain
        from anuga.file.sww import SWW_file
        from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import \
            Dirichlet_boundary

        Bd = Dirichlet_boundary([0.5, 0., 0.])

        # Create shallow water domain
        domain = Domain(*rectangular_cross(2, 2))
        domain.set_name('test1')
        domain.set_quantity('elevation', 2)
        domain.set_quantity('stage', 5)
        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
        for t in domain.evolve(yieldstep=0.5, finaltime=1):
            pass
            
        domain = Domain(*rectangular(3, 3))
        domain.set_name('test2')
        domain.set_quantity('elevation', 3)
        domain.set_quantity('stage', 50)
        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
        for t in domain.evolve(yieldstep=0.5, finaltime=1):
            pass
                
        outfile = 'test_out.sww'
        _sww_merge(['test1.sww', 'test2.sww'], outfile)
        self.assertTrue(os.access(outfile, os.F_OK))  
        
        # remove temp files
        if not sys.platform == 'win32':		
			os.remove('test1.sww')
			os.remove('test2.sww')
			os.remove(outfile)      

# 需要导入模块: from anuga.shallow_water.shallow_water_domain import Domain [as 别名]
# 或者: from anuga.shallow_water.shallow_water_domain.Domain import evolve [as 别名]
                            #width=25.0,
                            #depth=10.0,
                            #apron=5.0,
                            #manning=0.013,
                            #verbose=False)



#------------------------------------------------------------------------------
# Setup boundary conditions
#------------------------------------------------------------------------------
print 'Setting Boundary Conditions'
Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
Bi = anuga.Dirichlet_boundary([0.0, 0.0, 0.0])          # Inflow based on Flow Depth and Approaching Momentum !!!

Btus = anuga.Dirichlet_boundary([20.0, 0, 0])           # Outflow water at 10.0
Btds = anuga.Dirichlet_boundary([19.0, 0, 0])           # Outflow water at 9.0
domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br})



#------------------------------------------------------------------------------
# Evolve system through time
#------------------------------------------------------------------------------

for t in domain.evolve(yieldstep = 1, finaltime = 100):
    print domain.timestepping_statistics()

    domain.print_operator_timestepping_statistics()

# 需要导入模块: from anuga.shallow_water.shallow_water_domain import Domain [as 别名]
# 或者: from anuga.shallow_water.shallow_water_domain.Domain import evolve [as 别名]
    def test_read_sww(self):
        """
        Save to an sww file and then read back the info.
        Here we store the info "uniquely"
        """

        # ---------------------------------------------------------------------
        # Import necessary modules
        # ---------------------------------------------------------------------
        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
        from anuga.shallow_water.shallow_water_domain import Domain
        from anuga import Reflective_boundary
        from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import Dirichlet_boundary, Time_boundary

        # ---------------------------------------------------------------------
        # Setup computational domain
        # ---------------------------------------------------------------------
        length = 8.0
        width = 4.0
        dx = dy = 2  # Resolution: Length of subdivisions on both axes

        inc = 0.05  # Elevation increment

        points, vertices, boundary = rectangular_cross(int(length / dx), int(width / dy), len1=length, len2=width)
        domain = Domain(points, vertices, boundary)
        domain.set_name("read_sww_test" + str(domain.processor))  # Output name
        domain.set_quantities_to_be_stored({"elevation": 2, "stage": 2, "xmomentum": 2, "ymomentum": 2, "friction": 1})

        domain.set_store_vertices_uniquely(True)

        # ---------------------------------------------------------------------
        # Setup initial conditions
        # ---------------------------------------------------------------------
        domain.set_quantity("elevation", 0.0)  # Flat bed initially
        domain.set_quantity("friction", 0.01)  # Constant friction
        domain.set_quantity("stage", 0.0)  # Dry initial condition

        # ------------------------------------------------------------------
        # Setup boundary conditions
        # ------------------------------------------------------------------
        Bi = Dirichlet_boundary([0.4, 0, 0])  # Inflow
        Br = Reflective_boundary(domain)  # Solid reflective wall
        Bo = Dirichlet_boundary([-5, 0, 0])  # Outflow

        domain.set_boundary({"left": Bi, "right": Bo, "top": Br, "bottom": Br})

        # -------------------------------------------------------------------
        # Evolve system through time
        # -------------------------------------------------------------------

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

        # Check that quantities have been stored correctly
        source = domain.get_name() + ".sww"

        # x = fid.variables['x'][:]
        # y = fid.variables['y'][:]
        # stage = fid.variables['stage'][:]
        # elevation = fid.variables['elevation'][:]
        # fid.close()

        # assert len(stage.shape) == 2
        # assert len(elevation.shape) == 2

        # M, N = stage.shape

        sww_file = sww.Read_sww(source)

        # print 'last frame number',sww_file.get_last_frame_number()

        assert num.allclose(sww_file.x, domain.get_vertex_coordinates()[:, 0])
        assert num.allclose(sww_file.y, domain.get_vertex_coordinates()[:, 1])

        assert num.allclose(sww_file.time, [0.0, 1.0, 2.0, 3.0, 4.0])

        M = domain.get_number_of_triangles()

        assert num.allclose(num.reshape(num.arange(3 * M), (M, 3)), sww_file.vertices)

        last_frame_number = sww_file.get_last_frame_number()
        assert last_frame_number == 4

        assert num.allclose(sww_file.get_bounds(), [0.0, length, 0.0, width])

        assert "stage" in sww_file.quantities.keys()
        assert "friction" in sww_file.quantities.keys()
        assert "elevation" in sww_file.quantities.keys()
        assert "xmomentum" in sww_file.quantities.keys()
        assert "ymomentum" in sww_file.quantities.keys()

        for qname, q in sww_file.read_quantities(last_frame_number).items():

            # print qname
            # print num.linalg.norm(num.abs((domain.get_quantity(qname).get_values()-q).flatten()), ord=1)

            assert num.allclose(domain.get_quantity(qname).get_values(), q)

        # -----------------------------------------
        # Start the evolution off again at frame 3
#.........这里部分代码省略.........

# 需要导入模块: from anuga.shallow_water.shallow_water_domain import Domain [as 别名]
# 或者: from anuga.shallow_water.shallow_water_domain.Domain import evolve [as 别名]
    def test_get_mesh_and_quantities_from_unique_vertices_DE0_sww_file(self):
        """test_get_mesh_and_quantities_from_unique_vertices_sww_file(self):
        """     
        
        # Generate a test sww file with non trivial georeference
        
        import time, os

        # Setup
        #from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular

        # Create basic mesh (100m x 5m)
        width = 5
        length = 50
        t_end = 10
        points, vertices, boundary = rectangular(10, 1, length, width)

        # Create shallow water domain
        domain = Domain(points, vertices, boundary,
                        geo_reference = Geo_reference(56,308500,6189000))

        domain.set_name('test_get_mesh_and_quantities_from_unique_vertices_sww_file')
        swwfile = domain.get_name() + '.sww'
        domain.set_datadir('.')
        domain.set_flow_algorithm('DE0')
        domain.set_store_vertices_uniquely()

        Br = Reflective_boundary(domain)    # Side walls
        Bd = Dirichlet_boundary([1, 0, 0])  # inflow

        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})

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

        
        # Read it

        # Get mesh and quantities from sww file
        X = get_mesh_and_quantities_from_file(swwfile,
                                              quantities=['elevation',
                                                          'stage',
                                                          'xmomentum',
                                                          'ymomentum'], 
                                              verbose=False)
        mesh, quantities, time = X 

        #print quantities
        #print time

        dhash = domain.get_nodes()[:,0]*10+domain.get_nodes()[:,1]
        mhash = mesh.nodes[:,0]*10+mesh.nodes[:,1]        


        #print 'd_nodes',len(dhash)
        #print 'm_nodes',len(mhash)
        di = num.argsort(dhash)
        mi = num.argsort(mhash)
        minv = num.argsort(mi)
        dinv = num.argsort(di)

        #print 'd_tri',len(domain.get_triangles())
        #print 'm_tri',len(mesh.triangles)
        
        # Check that mesh has been recovered
        # triangle order should be ok
        assert num.allclose(mesh.nodes[mi,:],domain.get_nodes()[di,:])
        assert num.alltrue(minv[mesh.triangles] == dinv[domain.get_triangles()])


        # Check that time has been recovered
        assert num.allclose(time, range(t_end+1))

        z=domain.get_quantity('elevation').get_values(location='vertices').flatten()
        

        
        assert num.allclose(quantities['elevation'], z)

        for q in ['stage', 'xmomentum', 'ymomentum']:
            # Get quantity at last timestep
            q_ref=domain.get_quantity(q).get_values(location='vertices').flatten()

            #print q,quantities[q]
            q_sww=quantities[q][-1,:]
            
            msg = 'Quantity %s failed to be recovered' %q
            assert num.allclose(q_ref, q_sww, atol=1.0e-6), msg

# 需要导入模块: from anuga.shallow_water.shallow_water_domain import Domain [as 别名]
# 或者: from anuga.shallow_water.shallow_water_domain.Domain import evolve [as 别名]
    def test_get_mesh_and_quantities_from_de0_sww_file(self):
        """test_get_mesh_and_quantities_from_sww_file(self):
        """     
        
        # Generate a test sww file with non trivial georeference
        
        import time, os

        # Setup
        #from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular

        # Create basic mesh (100m x 5m)
        width = 5
        length = 50
        t_end = 10
        points, vertices, boundary = rectangular(length, width, 50, 5)

        # Create shallow water domain
        domain = Domain(points, vertices, boundary,
                        geo_reference = Geo_reference(56,308500,6189000))

        domain.set_name('test_get_mesh_and_quantities_from_sww_file')
        swwfile = domain.get_name() + '.sww'
        domain.set_datadir('.')
        domain.set_flow_algorithm('DE0')

        Br = Reflective_boundary(domain)    # Side walls
        Bd = Dirichlet_boundary([1, 0, 0])  # inflow

        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})

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

        
        # Read it

        # Get mesh and quantities from sww file
        X = get_mesh_and_quantities_from_file(swwfile,
                                              quantities=['elevation',
                                                          'stage',
                                                          'xmomentum',
                                                          'ymomentum'], 
                                              verbose=False)
        mesh, quantities, time = X
        

        
        # Check that mesh has been recovered
        assert num.alltrue(mesh.triangles == domain.get_triangles())
        assert num.allclose(mesh.nodes, domain.get_nodes())

        # Check that time has been recovered
        assert num.allclose(time, range(t_end+1))

        # Check that quantities have been recovered
        # (sww files use single precision)
        z=domain.get_quantity('elevation').get_values(location='unique vertices')
        assert num.allclose(quantities['elevation'], z)

        for q in ['stage', 'xmomentum', 'ymomentum']:
            # Get quantity at last timestep
            q_ref=domain.get_quantity(q).get_values(location='unique vertices')

            #print q,quantities[q]
            q_sww=quantities[q][-1,:]
            
            msg = 'Quantity %s failed to be recovered' %q
            assert num.allclose(q_ref, q_sww, atol=1.0e-2), msg

# 需要导入模块: from anuga.shallow_water.shallow_water_domain import Domain [as 别名]
# 或者: from anuga.shallow_water.shallow_water_domain.Domain import evolve [as 别名]
    def test_sww2domain1(self):
        ################################################
        #Create a test domain, and evolve and save it.
        ################################################
        #from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular

        #Create basic mesh

        yiel=0.01
        points, vertices, boundary = rectangular(10,10)

        #print "=============== boundary rect ======================="
        #print boundary

        #Create shallow water domain
        domain = Domain(points, vertices, boundary)
        domain.geo_reference = Geo_reference(56,11,11)
        domain.smooth = False
        domain.store = True
        domain.set_name('bedslope')
        domain.default_order=2
        #Bed-slope and friction
        domain.set_quantity('elevation', lambda x,y: -x/3)
        domain.set_quantity('friction', 0.1)
        # Boundary conditions
        from math import sin, pi
        Br = Reflective_boundary(domain)
        Bt = Transmissive_boundary(domain)
        Bd = Dirichlet_boundary([0.2,0.,0.])
        Bw = Time_boundary(domain=domain,function=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])

        #domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})

        domain.quantities_to_be_stored['xmomentum'] = 2
        domain.quantities_to_be_stored['ymomentum'] = 2
        #Initial condition
        h = 0.05
        elevation = domain.quantities['elevation'].vertex_values
        domain.set_quantity('stage', elevation + h)

        domain.check_integrity()
        #Evolution
        #domain.tight_slope_limiters = 1
        for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
            #domain.write_time()
            pass

        #print boundary


        filename = domain.datadir + os.sep + domain.get_name() + '.sww'
        domain2 = load_sww_as_domain(filename, None, fail_if_NaN=False,
                                        verbose=self.verbose)

        # Unfortunately we loss the boundaries top, bottom, left and right,
        # they are now all lumped into "exterior"

        #print "=============== boundary domain2 ======================="
        #print domain2.boundary
        

        #print domain2.get_boundary_tags()
        
        #points, vertices, boundary = rectangular(15,15)
        #domain2.boundary = boundary
        ###################
        ##NOW TEST IT!!!
        ###################

        os.remove(filename)

        bits = ['vertex_coordinates']
        for quantity in ['stage']:
            bits.append('get_quantity("%s").get_integral()' % quantity)
            bits.append('get_quantity("%s").get_values()' % quantity)

        for bit in bits:
            #print 'testing that domain.'+bit+' has been restored'
            #print bit
            #print 'done'
            #print eval('domain.'+bit)
            #print eval('domain2.'+bit)
            assert num.allclose(eval('domain.'+bit),eval('domain2.'+bit))

        ######################################
        #Now evolve them both, just to be sure
        ######################################x
        from time import sleep

        final = .1
        domain.set_quantity('friction', 0.1)
        domain.store = False
        domain.set_boundary({'exterior': Bd, 'left' : Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})


        for t in domain.evolve(yieldstep = yiel, finaltime = final):
            #domain.write_time()
            pass

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

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

        Test of get_maximum_inundation_elevation()
        and get_maximum_inundation_location().
   
        This is based on test_get_maximum_inundation_3(self) but works with the
        stored results instead of with the internal data structure.

        This test uses the underlying get_maximum_inundation_data for tests
        """

        verbose = False
        from anuga.config import minimum_storable_height
        
        initial_runup_height = -0.4
        final_runup_height = -0.3
        filename = 'runup_test_2'

        #--------------------------------------------------------------
        # Setup computational domain
        #--------------------------------------------------------------
        N = 10
        points, vertices, boundary = rectangular_cross(N, N)
        domain = Domain(points, vertices, boundary)
        domain.set_name(filename)
        domain.set_maximum_allowed_speed(1.0)
        #domain.set_minimum_storable_height(1.0e-5)
        domain.set_store_vertices_uniquely()

        # FIXME: This works better with old limiters so far
        domain.tight_slope_limiters = 0

        #--------------------------------------------------------------
        # Setup initial conditions
        #--------------------------------------------------------------
        def topography(x, y):
            return -x/2                             # linear bed slope

        # Use function for elevation
        domain.set_quantity('elevation', topography)
        domain.set_quantity('friction', 0.)                # Zero friction
        # Constant negative initial stage
        domain.set_quantity('stage', initial_runup_height)

        #--------------------------------------------------------------
        # Setup boundary conditions
        #--------------------------------------------------------------
        Br = Reflective_boundary(domain)                       # Reflective wall
        Bd = Dirichlet_boundary([final_runup_height, 0, 0])    # Constant inflow

        # All reflective to begin with (still water)
        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})

        #--------------------------------------------------------------
        # Test initial inundation height
        #--------------------------------------------------------------
        indices = domain.get_wet_elements()
        z = domain.get_quantity('elevation').\
                get_values(location='centroids', indices=indices)
        assert num.alltrue(z < initial_runup_height)

        q_ref = domain.get_maximum_inundation_elevation(minimum_height=minimum_storable_height)
        # First order accuracy
        assert num.allclose(q_ref, initial_runup_height, rtol=1.0/N)

        #--------------------------------------------------------------
        # Let triangles adjust
        #--------------------------------------------------------------
        q_max = None
        for t in domain.evolve(yieldstep = 0.1, finaltime = 1.0):
            q = domain.get_maximum_inundation_elevation(minimum_height=minimum_storable_height)

            if verbose:
                domain.write_time()
                print q
                
            if q > q_max:
                q_max = q

        #--------------------------------------------------------------
        # Test inundation height again
        #--------------------------------------------------------------
        #q_ref = domain.get_maximum_inundation_elevation()
        q = get_maximum_inundation_elevation(filename+'.sww')
        msg = 'We got %f, should have been %f' % (q, q_max)
        assert num.allclose(q, q_max, rtol=2.0/N), msg

        msg = 'We got %f, should have been %f' % (q, initial_runup_height)
        assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg

        # Test error condition if time interval is out
        try:
            q = get_maximum_inundation_elevation(filename+'.sww',
                                                 time_interval=[2.0, 3.0])
        except ValueError:
            pass
        else:
            msg = 'should have caught wrong time interval'
            raise Exception, msg
#.........这里部分代码省略.........

# 需要导入模块: from anuga.shallow_water.shallow_water_domain import Domain [as 别名]
# 或者: from anuga.shallow_water.shallow_water_domain.Domain import evolve [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 evolve [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 evolve [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()

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

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

        Test that the total flow through a cross section can be
        correctly obtained at run-time from the ANUGA domain.

        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
        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
        """

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

        for t in domain.evolve(yieldstep=0.1, finaltime=0.5):
            # Shortcuts to quantites
            stage = domain.get_quantity('stage')
            xmomentum = domain.get_quantity('xmomentum')
            ymomentum = domain.get_quantity('ymomentum')

            # 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 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)
                Q = domain.get_flow_through_cross_section(cross_section,
                                                          verbose=False)

                assert num.allclose(Q, uh*width)

        import cPickle        
        cPickle.dump(domain, open('domain_pickle.pickle', 'w'))
        domain_restored = cPickle.load(open('domain_pickle.pickle'))

        
        for t in domain_restored.evolve(yieldstep=0.1, finaltime=1.0):
            # Shortcuts to quantites
            stage = domain_restored.get_quantity('stage')
            xmomentum = domain_restored.get_quantity('xmomentum')
            ymomentum = domain_restored.get_quantity('ymomentum')
            
            # 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)
#.........这里部分代码省略.........

# 需要导入模块: from anuga.shallow_water.shallow_water_domain import Domain [as 别名]
# 或者: from anuga.shallow_water.shallow_water_domain.Domain import evolve [as 别名]

#.........这里部分代码省略.........
        xmomentum = fid.variables['xmomentum'][:]
        ymomentum = fid.variables['ymomentum'][:]
        elevation = fid.variables['elevation'][:]       

        # Create beginnings of boundary polygon based on sts_boundary
        boundary_polygon = create_sts_boundary(base_name)
        
        os.remove(order_file)

        # Append the remaining part of the boundary polygon to be defined by
        # the user
        bounding_polygon_utm=[]
        for point in bounding_polygon:
            zone,easting,northing=redfearn(point[0],point[1])
            bounding_polygon_utm.append([easting,northing])

        boundary_polygon.append(bounding_polygon_utm[3])
        boundary_polygon.append(bounding_polygon_utm[4])

        #print 'boundary_polygon', boundary_polygon
        
        plot=False
        if plot:
            from pylab import plot,show,axis
            boundary_polygon=ensure_numeric(boundary_polygon)
            bounding_polygon_utm=ensure_numeric(bounding_polygon_utm)
            #plot(lat_long_points[:,0],lat_long_points[:,1],'o')
            plot(boundary_polygon[:,0], boundary_polygon[:,1])
            plot(bounding_polygon_utm[:,0],bounding_polygon_utm[:,1])
            show()

        assert num.allclose(bounding_polygon_utm,boundary_polygon)


        extent_res=1000000
        meshname = 'urs_test_mesh' + '.tsh'
        interior_regions=None
        boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]}
        
        # have to change boundary tags from last example because now bounding
        # polygon starts in different place.
        create_mesh_from_regions(boundary_polygon,
                                 boundary_tags=boundary_tags,
                                 maximum_triangle_area=extent_res,
                                 filename=meshname,
                                 interior_regions=interior_regions,
                                 verbose=False)
        
        domain_fbound = Domain(meshname)
        domain_fbound.set_quantity('stage', tide)
        Bf = File_boundary(sts_file+'.sts', 
                           domain_fbound, 
                           boundary_polygon=boundary_polygon)
        Br = Reflective_boundary(domain_fbound)

        domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br})
        finaltime=time_step*(time_step_count-1)
        yieldstep=time_step
        temp_fbound=num.zeros(int(finaltime/yieldstep)+1,num.float)
    
        for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep,
                                                   finaltime=finaltime, 
                                                   skip_initial_step=False)):
            temp_fbound[i]=domain_fbound.quantities['stage'].centroid_values[2]
    
        
        domain_time = Domain(meshname)
        domain_time.set_quantity('stage', tide)
        Br = Reflective_boundary(domain_time)
        Bw = Time_boundary(domain=domain_time,
                         function=lambda t: [num.sin(t)+tide,3.*(20.+num.sin(t)+tide),2.*(20.+num.sin(t)+tide)])
        domain_time.set_boundary({'ocean': Bw,'otherocean': Br})
        
        temp_time=num.zeros(int(finaltime/yieldstep)+1,num.float)
        
        domain_time.set_starttime(domain_fbound.get_starttime())
        
        for i, t in enumerate(domain_time.evolve(yieldstep=yieldstep,
                                                   finaltime=finaltime, 
                                                   skip_initial_step=False)):
            temp_time[i]=domain_time.quantities['stage'].centroid_values[2]
        
        assert num.allclose(temp_fbound, temp_time)                
        assert num.allclose(domain_fbound.quantities['stage'].vertex_values,
                            domain_time.quantities['stage'].vertex_values)
                        
        assert num.allclose(domain_fbound.quantities['xmomentum'].vertex_values,
                            domain_time.quantities['xmomentum'].vertex_values)                        
                        
        assert num.allclose(domain_fbound.quantities['ymomentum'].vertex_values,
                            domain_time.quantities['ymomentum'].vertex_values)                                                
        

        try:
            os.remove(sts_file+'.sts')
        except:
            # Windoze can't remove this file for some reason 
            pass
        
        os.remove(meshname)