Python Domain.time方法代碼示例

# 需要導入模塊: from anuga.shallow_water.shallow_water_domain import Domain [as 別名]
# 或者: from anuga.shallow_water.shallow_water_domain.Domain import time [as 別名]
    def setUp(self):
        # print "****set up****"
        # Create an sww file

        # Set up an sww that has a geo ref.
        # have it cover an area in Australia.  'gong maybe
        # Don't have many triangles though!

        # Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
        # Zone:   56
        # Easting:  222908.705  Northing: 6233785.284
        # Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
        # Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660

        # geo-ref
        # Zone:   56
        # Easting:  220000  Northing: 6230000

        # have  a big area covered.
        mesh_file = tempfile.mktemp(".tsh")
        points_lat_long = [[-33, 152], [-35, 152], [-35, 150], [-33, 150]]
        spat = Geospatial_data(data_points=points_lat_long, points_are_lats_longs=True)
        points_ab = spat.get_data_points(absolute=True)

        geo = Geo_reference(56, 400000, 6000000)
        spat.set_geo_reference(geo)
        m = Mesh()
        m.add_vertices(spat)
        m.auto_segment()
        m.generate_mesh(verbose=False)
        m.export_mesh_file(mesh_file)

        # Create shallow water domain
        domain = Domain(mesh_file)

        os.remove(mesh_file)

        domain.default_order = 2
        # Set some field values
        # domain.set_quantity('stage', 1.0)
        domain.set_quantity("elevation", -0.5)
        domain.set_quantity("friction", 0.03)

        ######################
        # Boundary conditions
        B = Transmissive_boundary(domain)
        domain.set_boundary({"exterior": B})

        ######################
        # Initial condition - with jumps
        bed = domain.quantities["elevation"].vertex_values
        stage = num.zeros(bed.shape, num.float)

        h = 0.3
        for i in range(stage.shape[0]):
            if i % 2 == 0:
                stage[i, :] = bed[i, :] + h
            else:
                stage[i, :] = bed[i, :]

        domain.set_quantity("stage", stage)
        domain.set_quantity("xmomentum", stage * 22.0)
        domain.set_quantity("ymomentum", stage * 55.0)

        domain.distribute_to_vertices_and_edges()

        self.domain = domain

        C = domain.get_vertex_coordinates()
        self.X = C[:, 0:6:2].copy()
        self.Y = C[:, 1:6:2].copy()

        self.F = bed

        # sww_file = tempfile.mktemp("")
        self.domain.set_name("tid_P0")
        self.domain.format = "sww"
        self.domain.smooth = True
        self.domain.reduction = mean

        sww = SWW_file(self.domain)
        sww.store_connectivity()
        sww.store_timestep()
        self.domain.time = 2.0
        sww.store_timestep()
        self.sww = sww  # so it can be deleted

        # Create another sww file
        mesh_file = tempfile.mktemp(".tsh")
        points_lat_long = [[-35, 152], [-36, 152], [-36, 150], [-35, 150]]
        spat = Geospatial_data(data_points=points_lat_long, points_are_lats_longs=True)
        points_ab = spat.get_data_points(absolute=True)

        geo = Geo_reference(56, 400000, 6000000)
        spat.set_geo_reference(geo)
        m = Mesh()
        m.add_vertices(spat)
        m.auto_segment()
        m.generate_mesh(verbose=False)
        m.export_mesh_file(mesh_file)
#.........這裏部分代碼省略.........

# 需要導入模塊: from anuga.shallow_water.shallow_water_domain import Domain [as 別名]
# 或者: from anuga.shallow_water.shallow_water_domain.Domain import time [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 time [as 別名]
    def test_inundation_damage_list(self):

        # create mesh
        mesh_file = tempfile.mktemp(".tsh")
        points = [[0.0, 0.0], [6.0, 0.0], [6.0, 6.0], [0.0, 6.0]]
        m = Mesh()
        m.add_vertices(points)
        m.auto_segment()
        m.generate_mesh(verbose=False)
        m.export_mesh_file(mesh_file)

        # Create shallow water domain
        domain = Domain(mesh_file)
        os.remove(mesh_file)

        domain.default_order = 2

        # Set some field values
        domain.set_quantity("elevation", elevation_function)
        domain.set_quantity("friction", 0.03)
        domain.set_quantity("xmomentum", 22.0)
        domain.set_quantity("ymomentum", 55.0)

        ######################
        # Boundary conditions
        B = Transmissive_boundary(domain)
        domain.set_boundary({"exterior": B})

        # This call mangles the stage values.
        domain.distribute_to_vertices_and_edges()
        domain.set_quantity("stage", 0.3)

        # sww_file = tempfile.mktemp("")
        domain.set_name("datatest" + str(time.time()))
        domain.format = "sww"
        domain.smooth = True
        domain.reduction = mean

        sww = SWW_file(domain)
        sww.store_connectivity()
        sww.store_timestep()
        domain.set_quantity("stage", -0.3)
        domain.time = 2.0
        sww.store_timestep()

        # Create a csv file
        csv_file = tempfile.mktemp(".csv")
        fd = open(csv_file, "wb")
        writer = csv.writer(fd)
        writer.writerow(["x", "y", STR_VALUE_LABEL, CONT_VALUE_LABEL, "ROOF_TYPE", WALL_TYPE_LABEL, SHORE_DIST_LABEL])
        writer.writerow([5.5, 0.5, "10", "130000", "Metal", "Timber", 20])
        writer.writerow([4.5, 1.0, "150", "76000", "Metal", "Double Brick", 20])
        writer.writerow([0.1, 1.5, "100", "76000", "Metal", "Brick Veneer", 300])
        writer.writerow([6.1, 1.5, "100", "76000", "Metal", "Brick Veneer", 300])
        fd.close()

        extension = ".csv"
        csv_fileII = tempfile.mktemp(extension)
        fd = open(csv_fileII, "wb")
        writer = csv.writer(fd)
        writer.writerow(["x", "y", STR_VALUE_LABEL, CONT_VALUE_LABEL, "ROOF_TYPE", WALL_TYPE_LABEL, SHORE_DIST_LABEL])
        writer.writerow([5.5, 0.5, "10", "130000", "Metal", "Timber", 20])
        writer.writerow([4.5, 1.0, "150", "76000", "Metal", "Double Brick", 20])
        writer.writerow([0.1, 1.5, "100", "76000", "Metal", "Brick Veneer", 300])
        writer.writerow([6.1, 1.5, "100", "76000", "Metal", "Brick Veneer", 300])
        fd.close()

        sww_file = domain.get_name() + "." + domain.format
        # print "sww_file",sww_file
        marker = "_gosh"
        inundation_damage(sww_file, [csv_file, csv_fileII], exposure_file_out_marker=marker, verbose=False)

        # Test one file
        csv_handle = Exposure(csv_file[:-4] + marker + extension)
        struct_loss = csv_handle.get_column(EventDamageModel.STRUCT_LOSS_TITLE)
        # print "struct_loss",struct_loss
        struct_loss = [float(x) for x in struct_loss]
        # pprint(struct_loss)
        assert num.allclose(struct_loss, [10.0, 150.0, 66.55333347876866, 0.0])
        depth = csv_handle.get_column(EventDamageModel.MAX_DEPTH_TITLE)
        # print "depth",depth
        depth = [float(x) for x in depth]
        assert num.allclose(depth, [3.000000011920929, 2.9166666785875957, 2.2666666785875957, -0.3])

        # Test another file
        csv_handle = Exposure(csv_fileII[:-4] + marker + extension)
        struct_loss = csv_handle.get_column(EventDamageModel.STRUCT_LOSS_TITLE)
        # print "struct_loss",struct_loss
        struct_loss = [float(x) for x in struct_loss]

        # pprint(struct_loss)
        assert num.allclose(struct_loss, [10.0, 150.0, 66.553333478768664, 0.0])
        depth = csv_handle.get_column(EventDamageModel.MAX_DEPTH_TITLE)
        # print "depth",depth
        depth = [float(x) for x in depth]
        assert num.allclose(depth, [3.000000011920929, 2.9166666785875957, 2.2666666785875957, -0.3])
        os.remove(sww.filename)
        os.remove(csv_file)
        os.remove(csv_fileII)