Python Domain.set_boundary方法代码示例

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import set_boundary [as 别名]
    def test_rate_operator_negative_rate_full(self):
        from anuga.config import rho_a, rho_w, eta_w
        from math import pi, cos, sin

        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)

        # Flat surface with 1m of water
        domain.set_quantity("elevation", 0)
        domain.set_quantity("stage", 10.0)
        domain.set_quantity("friction", 0)

        Br = Reflective_boundary(domain)
        domain.set_boundary({"exterior": Br})

        #        print domain.quantities['elevation'].centroid_values
        #        print domain.quantities['stage'].centroid_values
        #        print domain.quantities['xmomentum'].centroid_values
        #        print domain.quantities['ymomentum'].centroid_values

        # Apply operator to these triangles
        indices = [0, 1, 3]

        # Catchment_Rain_Polygon = read_polygon(join('CatchmentBdy.csv'))
        # rainfall = file_function(join('1y120m.tms'), quantities=['rainfall'])
        rate = -1.0
        factor = 10.0
        default_rate = 0.0

        operator = Rate_operator(domain, rate=rate, factor=factor, indices=None, default_rate=default_rate)

        # Apply Operator
        domain.timestep = 2.0
        operator()

        stage_ex = [0.0, 0.0, 0.0, 0.0]
        step_integral = -80.0

        # print domain.quantities['elevation'].centroid_values
        # print domain.quantities['stage'].centroid_values
        # print domain.quantities['xmomentum'].centroid_values
        # print domain.quantities['ymomentum'].centroid_values
        # print domain.fractional_step_volume_integral

        assert num.allclose(domain.quantities["stage"].centroid_values, stage_ex)
        assert num.allclose(domain.quantities["xmomentum"].centroid_values, 0.0)
        assert num.allclose(domain.quantities["ymomentum"].centroid_values, 0.0)
        assert num.allclose(domain.fractional_step_volume_integral, step_integral)

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import set_boundary [as 别名]
    def test_set_w_uh_vh_operator_time(self):
        from anuga.config import rho_a, rho_w, eta_w
        from math import pi, cos, sin

        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)

        #Flat surface with 1m of water
        domain.set_quantity('elevation', 0)
        domain.set_quantity('stage', 1.0)
        domain.set_quantity('xmomentum', 7.0)
        domain.set_quantity('ymomentum', 8.0)
        domain.set_quantity('friction', 0)

        Br = Reflective_boundary(domain)
        domain.set_boundary({'exterior': Br})


#        print domain.quantities['w_uh_vh'].centroid_values
#        print domain.quantities['xmomentum'].centroid_values
#        print domain.quantities['ymomentum'].centroid_values

        # Apply operator to these triangles
        indices = [0,1,3]

        w_uh_vh = lambda t : [t, t+1, t+2]


        operator = Set_w_uh_vh_operator(domain, w_uh_vh=w_uh_vh, indices=indices)
        
        # Apply Operator
        domain.timestep = 2.0
        domain.time = 1.0
        operator()

        t = domain.time
        stage_ex = [ t,  t,   1.,  t]
        xmom_ex = [ t+1,  t+1,   7.,  t+1]
        ymom_ex = [ t+2,  t+2,   8.,  t+2]


        #print domain.quantities['stage'].centroid_values
        #print domain.quantities['xmomentum'].centroid_values
        #print domain.quantities['ymomentum'].centroid_values

        assert num.allclose(domain.quantities['stage'].centroid_values, stage_ex)
        assert num.allclose(domain.quantities['xmomentum'].centroid_values, xmom_ex)
        assert num.allclose(domain.quantities['ymomentum'].centroid_values, ymom_ex)

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import set_boundary [as 别名]
    def test_set_stage_operator_negative(self):
        from anuga.config import rho_a, rho_w, eta_w
        from math import pi, cos, sin

        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)

        #Flat surface with 1m of water
        domain.set_quantity('elevation', lambda x,y : -2*x)
        domain.set_quantity('stage', 1.0)
        domain.set_quantity('friction', 0)

        Br = Reflective_boundary(domain)
        domain.set_boundary({'exterior': Br})

#        print domain.quantities['elevation'].centroid_values
#        print domain.quantities['stage'].centroid_values
#        print domain.quantities['xmomentum'].centroid_values
#        print domain.quantities['ymomentum'].centroid_values

        # Apply operator to these triangles
        indices = [0,1,3]



        #Catchment_Rain_Polygon = read_polygon(join('CatchmentBdy.csv'))
        #rainfall = file_function(join('1y120m.tms'), quantities=['rainfall'])
        stage = -5.0


        operator = Set_stage_operator(domain, stage=stage, indices=indices)


        # Apply Operator
        domain.timestep = 2.0
        operator()

        stage_ex = [ -5.,  -5.,   1.,  -5.]

        #print domain.quantities['elevation'].centroid_values
        #print domain.quantities['stage'].centroid_values
        #print domain.quantities['xmomentum'].centroid_values
        #print domain.quantities['ymomentum'].centroid_values

        assert num.allclose(domain.quantities['stage'].centroid_values, stage_ex)
        assert num.allclose(domain.quantities['xmomentum'].centroid_values, 0.0)
        assert num.allclose(domain.quantities['ymomentum'].centroid_values, 0.0)

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import set_boundary [as 别名]
    def test_rate_operator_simple(self):
        from anuga.config import rho_a, rho_w, eta_w
        from math import pi, cos, sin

        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)

        #Flat surface with 1m of water
        domain.set_quantity('elevation', 0)
        domain.set_quantity('stage', 1.0)
        domain.set_quantity('friction', 0)

        Br = Reflective_boundary(domain)
        domain.set_boundary({'exterior': Br})


#        print domain.quantities['stage'].centroid_values
#        print domain.quantities['xmomentum'].centroid_values
#        print domain.quantities['ymomentum'].centroid_values

        # Apply operator to these triangles
        indices = [0,1,3]

        rate = 1.0
        factor = 10.0
        default_rate= 0.0

        operator = Rate_operator(domain, rate=rate, factor=factor, \
                      indices=indices, default_rate = default_rate)
        
        # Apply Operator
        domain.timestep = 2.0
        operator()

        stage_ex = [ 21.,  21.,   1.,  21.]

#        print domain.quantities['stage'].centroid_values
#        print domain.quantities['xmomentum'].centroid_values
#        print domain.quantities['ymomentum'].centroid_values

        assert num.allclose(domain.quantities['stage'].centroid_values, stage_ex)
        assert num.allclose(domain.quantities['xmomentum'].centroid_values, 0.0)
        assert num.allclose(domain.quantities['ymomentum'].centroid_values, 0.0)
        assert num.allclose(domain.fractional_step_volume_integral, factor*domain.timestep*(rate*domain.areas[indices]).sum())

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import set_boundary [as 别名]
    def test_set_stage_operator_simple(self):
        from anuga.config import rho_a, rho_w, eta_w
        from math import pi, cos, sin

        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)

        #Flat surface with 1m of water
        domain.set_quantity('elevation', 0)
        domain.set_quantity('stage', 1.0)
        domain.set_quantity('friction', 0)

        Br = Reflective_boundary(domain)
        domain.set_boundary({'exterior': Br})


#        print domain.quantities['stage'].centroid_values
#        print domain.quantities['xmomentum'].centroid_values
#        print domain.quantities['ymomentum'].centroid_values

        # Apply operator to these triangles
        indices = [0,1,3]

        stage = 3.0


        operator = Set_stage_operator(domain, stage=stage, indices=indices)
        
        # Apply Operator
        domain.timestep = 2.0
        operator()

        stage_ex = [ 3.,  3.,   1.,  3.]


        #print domain.quantities['stage'].centroid_values
        #print domain.quantities['xmomentum'].centroid_values
        #print domain.quantities['ymomentum'].centroid_values

        assert num.allclose(domain.quantities['stage'].centroid_values, stage_ex)
        assert num.allclose(domain.quantities['xmomentum'].centroid_values, 0.0)
        assert num.allclose(domain.quantities['ymomentum'].centroid_values, 0.0)

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import set_boundary [as 别名]
    def concept_ungenerateIII(self):
        from anuga import Domain, Reflective_boundary, \
                            Dirichlet_boundary
        from anuga.pmesh.mesh_interface import create_mesh_from_regions

        # These are the absolute values
        polygon = [[0,0], [100,0], [100,100], [0,100]]

        boundary_tags = {'wall': [0,1,3], 'wave': [2]}
        inner1_polygon = [[10,10], [20,10], [20,20], [10,20]]
        inner2_polygon = [[30,30], [40,30], [40,40], [30,40]]

        max_area = 1
        interior_regions = [(inner1_polygon, 5), (inner2_polygon, 10)]
        m = create_mesh_from_regions(polygon,
                                     boundary_tags,
                                     max_area,
                                     interior_regions=interior_regions)

        fileName = tempfile.mktemp('.txt')
        file = open(fileName, 'w')
        file.write('         1       ??      ??\n\
       90.0       90.0\n\
       81.0       90.0\n\
       81.0       81.0\n\
       90.0       81.0\n\
       90.0       90.0\n\
END\n\
         2      ?? ??\n\
       10.0       80.0\n\
       10.0       90.0\n\
       20.0       90.0\n\
       10.0       80.0\n\
END\n\
END\n')
        file.close()

        m.import_ungenerate_file(fileName)
        os.remove(fileName)
        m.generate_mesh(maximum_triangle_area=max_area, verbose=False)
        mesh_filename = 'mesh.tsh'
        m.export_mesh_file(mesh_filename)

        domain = Domain(mesh_filename, use_cache=False)

        Br = Reflective_boundary(domain)
        Bd = Dirichlet_boundary([3, 0, 0])
        domain.set_boundary({'wall': Br, 'wave': Bd})
        yieldstep = 0.1
        finaltime = 10
        for t in domain.evolve(yieldstep, finaltime):
            domain.write_time()

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import set_boundary [as 别名]
    def test_rate_operator_functions_rate_default_rate(self):
        from anuga.config import rho_a, rho_w, eta_w
        from math import pi, cos, sin

        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)

        #Flat surface with 1m of water
        domain.set_quantity('elevation', 0)
        domain.set_quantity('stage', 1.0)
        domain.set_quantity('friction', 0)

        Br = Reflective_boundary(domain)
        domain.set_boundary({'exterior': Br})

        verbose = False
        
        if verbose:
            print domain.quantities['elevation'].centroid_values
            print domain.quantities['stage'].centroid_values
            print domain.quantities['xmomentum'].centroid_values
            print domain.quantities['ymomentum'].centroid_values

        # Apply operator to these triangles
        indices = [0,1,3]
        factor = 10.0


        def main_rate(t):
            if t > 20:
                msg = 'Model time exceeded.'
                raise Modeltime_too_late, msg
            else:
                return 3.0 * t + 7.0

        default_rate = lambda t: 3*t + 7


        operator = Rate_operator(domain, rate=main_rate, factor=factor, \
                      indices=indices, default_rate = default_rate)


        # Apply Operator
        domain.timestep = 2.0
        operator()

        t = operator.get_time()
        d = operator.get_timestep()*main_rate(t)*factor + 1
        stage_ex = [ d,  d,   1.,  d]

        if verbose:
            print domain.quantities['elevation'].centroid_values
            print domain.quantities['stage'].centroid_values
            print domain.quantities['xmomentum'].centroid_values
            print domain.quantities['ymomentum'].centroid_values

        assert num.allclose(domain.quantities['stage'].centroid_values, stage_ex)
        assert num.allclose(domain.quantities['xmomentum'].centroid_values, 0.0)
        assert num.allclose(domain.quantities['ymomentum'].centroid_values, 0.0)
        assert num.allclose(domain.fractional_step_volume_integral, ((d-1.)*domain.areas[indices]).sum())

        domain.set_starttime(30.0)
        domain.timestep = 1.0
        operator()

        t = operator.get_time()
        d = operator.get_timestep()*default_rate(t)*factor + d
        stage_ex = [ d,  d,   1.,  d]

        if verbose:
            print domain.quantities['elevation'].centroid_values
            print domain.quantities['stage'].centroid_values
            print domain.quantities['xmomentum'].centroid_values
            print domain.quantities['ymomentum'].centroid_values

        assert num.allclose(domain.quantities['stage'].centroid_values, stage_ex)
        assert num.allclose(domain.quantities['xmomentum'].centroid_values, 0.0)
        assert num.allclose(domain.quantities['ymomentum'].centroid_values, 0.0)

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import set_boundary [as 别名]
    def test_rate_operator_rate_from_file(self):
        from anuga.config import rho_a, rho_w, eta_w
        from math import pi, cos, sin

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


        #---------------------------------
        #Typical ASCII file
        #---------------------------------
        finaltime = 1200
        filename = 'test_file_function'
        fid = open(filename + '.txt', 'w')
        start = time.mktime(time.strptime('2000', '%Y'))
        dt = 60  #One minute intervals
        t = 0.0
        while t <= finaltime:
            t_string = time.strftime(time_format, time.gmtime(t+start))
            fid.write('%s, %f %f %f\n' %(t_string, 2*t, t**2, sin(t*pi/600)))
            t += dt

        fid.close()

        #Convert ASCII file to NetCDF (Which is what we really like!)
        timefile2netcdf(filename+'.txt')


        #Create file function from time series
        F = file_function(filename + '.tms',
                          quantities = ['Attribute0',
                                        'Attribute1',
                                        'Attribute2'])


        #Now try interpolation
        for i in range(20):
            t = i*10
            q = F(t)

            #Exact linear intpolation
            assert num.allclose(q[0], 2*t)
            if i%6 == 0:
                assert num.allclose(q[1], t**2)
                assert num.allclose(q[2], sin(t*pi/600))

        #Check non-exact

        t = 90 #Halfway between 60 and 120
        q = F(t)
        assert num.allclose( (120**2 + 60**2)/2, q[1] )
        assert num.allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] )


        t = 100 #Two thirds of the way between between 60 and 120
        q = F(t)
        assert num.allclose( 2*120**2/3 + 60**2/3, q[1] )
        assert num.allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )

        #os.remove(filename + '.txt')
        #os.remove(filename + '.tms')


        domain = Domain(points, vertices)

        #Flat surface with 1m of water
        domain.set_quantity('elevation', 0)
        domain.set_quantity('stage', 1.0)
        domain.set_quantity('friction', 0)

        Br = Reflective_boundary(domain)
        domain.set_boundary({'exterior': Br})

#        print domain.quantities['elevation'].centroid_values
#        print domain.quantities['stage'].centroid_values
#        print domain.quantities['xmomentum'].centroid_values
#        print domain.quantities['ymomentum'].centroid_values

        # Apply operator to these triangles
        indices = [0,1,3]


        rate = file_function(filename + '.tms', quantities=['Attribute1'])
        

        # Make starttime of domain consistent with tms file starttime
        domain.set_starttime(rate.starttime)
                    
        factor = 1000.0
        default_rate= 17.7

        operator = Rate_operator(domain, rate=rate, factor=factor, \
#.........这里部分代码省略.........

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import set_boundary [as 别名]
    def concept_ungenerateII(self):
        from anuga import Domain, Reflective_boundary, Dirichlet_boundary

        x=0
        y=0
        mesh_geo = geo_reference=Geo_reference(56, x, y)

        # These are the absolute values
        polygon_absolute = [[0,0], [100,0], [100,100], [0,100]]
        x_p = -10
        y_p = -40
        geo_ref_poly = Geo_reference(56, x_p, y_p)
        polygon = geo_ref_poly.change_points_geo_ref(polygon_absolute)

        boundary_tags = {'wall': [0,1,3], 'wave': [2]}

        inner1_polygon_absolute = [[10,10], [20,10], [20,20], [10,20]]
        inner1_polygon = geo_ref_poly.\
                            change_points_geo_ref(inner1_polygon_absolute)

        inner2_polygon_absolute = [[30,30], [40,30], [40,40], [30,40]]
        inner2_polygon = geo_ref_poly.\
                            change_points_geo_ref(inner2_polygon_absolute)

        max_area = 1
        interior_regions = [(inner1_polygon, 5), (inner2_polygon, 10)]
        m = create_mesh_from_regions(polygon,
                                     boundary_tags,
                                     max_area,
                                     interior_regions=interior_regions,
                                     poly_geo_reference=geo_ref_poly,
                                     mesh_geo_reference=mesh_geo)

        m.export_mesh_file('a_test_mesh_iknterface.tsh')

        fileName = tempfile.mktemp('.txt')
        file = open(fileName, 'w')
        file.write('         1       ??      ??\n\
       90.0       90.0\n\
       81.0       90.0\n\
       81.0       81.0\n\
       90.0       81.0\n\
       90.0       90.0\n\
END\n\
         2      ?? ??\n\
       10.0       80.0\n\
       10.0       90.0\n\
       20.0       90.0\n\
       10.0       80.0\n\
END\n\
END\n')
        file.close()

        m.import_ungenerate_file(fileName)      #, tag='wall')
        os.remove(fileName)
        m.generate_mesh(maximum_triangle_area=max_area, verbose=False)
        mesh_filename = 'bento_b.tsh'
        m.export_mesh_file(mesh_filename)

        domain = Domain(mesh_filename, use_cache = False)

        Br = Reflective_boundary(domain)
        Bd = Dirichlet_boundary([3, 0, 0])
        domain.set_boundary({'wall': Br, 'wave': Bd})
        yieldstep = 0.1
        finaltime = 10
        for t in domain.evolve(yieldstep, finaltime):
            domain.write_time()

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import set_boundary [as 别名]
    def test_set_quantity_negative(self):
        from anuga.config import rho_a, rho_w, eta_w
        from math import pi, cos, sin

        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)

        # Flat surface with 1m of water
        domain.set_quantity("elevation", lambda x, y: -2 * x)
        domain.set_quantity("stage", 1.0)
        domain.set_quantity("friction", 0)

        Br = Reflective_boundary(domain)
        domain.set_boundary({"exterior": Br})

        #        print domain.quantities['elevation'].centroid_values
        #        print domain.quantities['stage'].centroid_values
        #        print domain.quantities['xmomentum'].centroid_values
        #        print domain.quantities['ymomentum'].centroid_values

        # Apply operator to these triangles
        indices = [0, 1, 3]

        # Catchment_Rain_Polygon = read_polygon(join('CatchmentBdy.csv'))
        # rainfall = file_function(join('1y120m.tms'), quantities=['rainfall'])
        stage = -5.0

        try:
            update_stage = Set_quantity(domain, "stage", value=stage, indices=indices)
        except AssertionError:
            pass
        except e:
            self.fail("Unexpected exception thrown:", e)
        else:
            self.fail("ExpectedException not thrown")

        update_stage = Set_quantity(domain, "stage", value=stage, indices=indices, test_stage=False)
        # Apply Operator
        domain.timestep = 2.0
        update_stage()

        stage_ex = [-5.0, -5.0, 1.0, -5.0]

        # print domain.quantities['elevation'].centroid_values
        # print domain.quantities['stage'].centroid_values
        # print domain.quantities['xmomentum'].centroid_values
        # print domain.quantities['ymomentum'].centroid_values

        assert num.allclose(domain.quantities["stage"].centroid_values, stage_ex)
        assert num.allclose(domain.quantities["xmomentum"].centroid_values, 0.0)
        assert num.allclose(domain.quantities["ymomentum"].centroid_values, 0.0)

        update_stage = Set_stage(domain, stage=stage, indices=indices)
        domain.timestep = 2.0
        update_stage()

        stage_ex = [-1.33333333, -2.66666667, 1.0, -1.33333333]

        assert num.allclose(domain.quantities["stage"].centroid_values, stage_ex)
        assert num.allclose(domain.quantities["xmomentum"].centroid_values, 0.0)
        assert num.allclose(domain.quantities["ymomentum"].centroid_values, 0.0)

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import set_boundary [as 别名]
    def test_runup_sinusoid(self):
        """ Run a version of the validation test runup_sinusoid
        to ensure limiting solution has small velocity
        """

        points, vertices, boundary = anuga.rectangular_cross(20,20, len1=1., len2=1.)


        domain=Domain(points,vertices,boundary)    # Create Domain
        domain.set_flow_algorithm('DE0')
        
        domain.set_name('runup_sinusoid_v2')                         # Output to file runup.sww
        domain.set_datadir('.')                          # Use current folder
        domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
        #domain.set_store_vertices_uniquely(True)
        #------------------
        # Define topography
        #------------------
        scale_me=1.0

        def topography(x,y):
            return (-x/2.0 +0.05*num.sin((x+y)*50.0))*scale_me

        def stagefun(x,y):
            stge=-0.2*scale_me #+0.01*(x>0.9)
            return stge

        domain.set_quantity('elevation',topography)     # Use function for elevation
        domain.get_quantity('elevation').smooth_vertex_values()
        domain.set_quantity('friction',0.03)            # Constant friction


        domain.set_quantity('stage', stagefun)             # Constant negative initial stage
        domain.get_quantity('stage').smooth_vertex_values()


        #--------------------------
        # Setup boundary conditions
        #--------------------------
        Br=anuga.Reflective_boundary(domain)                 # Solid reflective wall
        Bd=anuga.Dirichlet_boundary([-0.1*scale_me,0.,0.])   # Constant boundary values -- not used in this example

        #----------------------------------------------
        # Associate boundary tags with boundary objects
        #----------------------------------------------
        domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom':Br})

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

        for t in domain.evolve(yieldstep=7.0,finaltime=7.0):
            #print domain.timestepping_statistics()
            xx = domain.quantities['xmomentum'].centroid_values
            yy = domain.quantities['ymomentum'].centroid_values
            dd = domain.quantities['stage'].centroid_values - domain.quantities['elevation'].centroid_values
            #dd_raw=1.0*dd
            dd = (dd)*(dd>1.0e-03)+1.0e-03
            vv = ( (xx/dd)**2 + (yy/dd)**2)**0.5
            vv = vv*(dd>1.0e-03)
            #print 'Peak velocity is: ', vv.max(), vv.argmax()
            #print 'Volume is', sum(dd_raw*domain.areas)


        #print vv.max()

        assert num.all(vv<1.01e-01)

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import set_boundary [as 别名]
# -----------------------------------------------------------------------------
domain.forcing_terms.append(Coulomb_forcing_term)


# -----------------------------------------------------------------------------
# Setup boundary conditions
# ------------------------------------------------------------------------------
from math import sin, pi, exp

Br = anuga.Reflective_boundary(domain)  # Solid reflective wall
Bt = anuga.Transmissive_boundary(domain)  # Continue all values on boundary
# Bd = anuga.Dirichlet_boundary([1,0.,0.]) # Constant boundary values
BTime = anuga.Time_boundary(domain, f_right)

# Associate boundary tags with boundary objects
domain.set_boundary({"left": Bt, "right": BTime, "top": Br, "bottom": Br})


# ------------------------------------------------------------------------------
# Produce a documentation of parameters
# ------------------------------------------------------------------------------
if myid == 0:
    parameter_file = open("parameters.tex", "w")
    parameter_file.write("\\begin{verbatim}\n")
    from pprint import pprint

    pprint(domain.get_algorithm_parameters(), parameter_file, indent=4)
    parameter_file.write("\\end{verbatim}\n")
    parameter_file.close()

# ------------------------------------------------------------------------------

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import set_boundary [as 别名]
    def test_rate_operator_rate_quantity(self):
        from anuga.config import rho_a, rho_w, eta_w
        from math import pi, cos, sin

        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)

        #Flat surface with 1m of water
        domain.set_quantity('elevation', 0.0)
        domain.set_quantity('stage', 1.0)
        domain.set_quantity('friction', 0.0)

        Br = Reflective_boundary(domain)
        domain.set_boundary({'exterior': Br})

        verbose = False

        if verbose:
            print domain.quantities['elevation'].centroid_values
            print domain.quantities['stage'].centroid_values
            print domain.quantities['xmomentum'].centroid_values
            print domain.quantities['ymomentum'].centroid_values

        # Apply operator to these triangles
        indices = [0,1,3]
        factor = 10.0


        from anuga import Quantity
        rate_Q = Quantity(domain)
        rate_Q.set_values(1.0)

        operator = Rate_operator(domain, rate=rate_Q, factor=factor, \
                                 indices=indices)


        # Apply Operator
        domain.timestep = 2.0
        operator()
        rate = rate_Q.centroid_values[indices]
        t = operator.get_time()
        Q = operator.get_Q()

        rate = rate*factor
        Q_ex = num.sum(domain.areas[indices]*rate)
        d = operator.get_timestep()*rate + 1


        #print "d"
        #print d
        #print Q_ex
        #print Q
        stage_ex = num.array([ 1.0,  1.0,   1.0,  1.0])
        stage_ex[indices] = d
        
        verbose = False
        
        if verbose:
            print domain.quantities['elevation'].centroid_values
            print domain.quantities['stage'].centroid_values
            print domain.quantities['xmomentum'].centroid_values
            print domain.quantities['ymomentum'].centroid_values

        assert num.allclose(domain.quantities['stage'].centroid_values, stage_ex)
        assert num.allclose(domain.quantities['xmomentum'].centroid_values, 0.0)
        assert num.allclose(domain.quantities['ymomentum'].centroid_values, 0.0)
        assert num.allclose(Q_ex, Q)
        assert num.allclose(domain.fractional_step_volume_integral, ((d-1.)*domain.areas[indices]).sum())

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import set_boundary [as 别名]
    def test_rate_operator_functions_spatial(self):
        from anuga.config import rho_a, rho_w, eta_w
        from math import pi, cos, sin

        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)

        area = numpy.sum(domain.areas)

        # Flat surface with 1m of water
        domain.set_quantity("elevation", 0.0)
        domain.set_quantity("stage", 1.0)
        domain.set_quantity("friction", 0.0)

        Br = Reflective_boundary(domain)
        domain.set_boundary({"exterior": Br})

        verbose = False

        if verbose:
            print domain.quantities["elevation"].centroid_values
            print domain.quantities["stage"].centroid_values
            print domain.quantities["xmomentum"].centroid_values
            print domain.quantities["ymomentum"].centroid_values

        # Apply operator to these triangles
        factor = 10.0

        def main_spatial_rate(x, y, t):
            # x and y should be an n by 1 array
            return x + y

        default_rate = 0.0

        operator = Rate_operator(domain, rate=main_spatial_rate, factor=factor, default_rate=default_rate)

        # Apply Operator
        domain.timestep = 2.0
        operator()

        t = operator.get_time()
        Q = operator.get_Q()
        x = operator.coord_c[:, 0]
        y = operator.coord_c[:, 1]
        rate = main_spatial_rate(x, y, t) * factor
        Q_ex = num.sum(domain.areas * rate)
        d = operator.get_timestep() * rate + 1

        # print "d"
        # print d
        # print area, Q, Q_ex
        stage_ex = num.array([1.0, 1.0, 1.0, 1.0])
        stage_ex[:] = d

        if verbose:
            print domain.quantities["elevation"].centroid_values
            print domain.quantities["stage"].centroid_values
            print domain.quantities["xmomentum"].centroid_values
            print domain.quantities["ymomentum"].centroid_values

        assert num.allclose(domain.quantities["stage"].centroid_values, stage_ex)
        assert num.allclose(domain.quantities["xmomentum"].centroid_values, 0.0)
        assert num.allclose(domain.quantities["ymomentum"].centroid_values, 0.0)
        assert num.allclose(Q_ex, Q)

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

domain.set_quantity("concentration", 0.01)
domain.set_quantity("elevation", topography)  # elevation is a function
domain.set_quantity("friction", 0.01)  # Constant friction
domain.set_quantity("stage", topography)  # Dry initial condition

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

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

# ------------------------------------------------------------------------------
# Setup erosion operator in the middle of dam
# ------------------------------------------------------------------------------
print "Set up Erosion Area to test..."

from anuga.operators.sed_transport_operator import Sed_transport_operator

# create operator
op1 = Sed_transport_operator(domain)
# op1.set_inflow_concentration(0.02)


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