Python Domain.set_quantity方法代码示例

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import set_quantity [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_quantity [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_quantity [as 别名]
    def test_slide_tsunami_domain(self):

        if anuga_installed:
            pass
        else:
            print 'Note: test_slide_tsunami_domain not tested as ANUGA '\
                'is not installed'
            return
        length = 600.0
        dep = 150.0
        th = 9.0
        thk = 15.0
        wid = 340.0
        kappa = 3.0
        kappad = 0.8
        x0 = 100000.
        y0 = x0
        
        from anuga.pmesh.mesh_interface import create_mesh_from_regions
        polygon = [[0,0],[200000,0],[200000,200000],[0,200000]]
        create_mesh_from_regions(polygon,
                                 {'e0': [0], 'e1': [1], 'e2': [2], 'e3': [3]},
                                 maximum_triangle_area=5000000000,
                                 filename='test.msh',
                                 verbose = False)

        domain = Domain('test.msh', use_cache = True, verbose = False)

        slide = slide_tsunami(length, dep, th, x0, y0, \
                              wid, thk, kappa, kappad, \
                              domain=domain,verbose=False)

        domain.set_quantity('stage', slide)
        stage = domain.get_quantity('stage')
        w = stage.get_values()

##        check = [[-0.0 -0.0 -0.0],
##                 [-.189709745 -517.877716 -0.0],
##                 [-0.0 -0.0 -2.7695931e-08],
##                 [-0.0 -2.7695931e-08 -1.897097e-01]
##                 [-0.0 -517.877716 -0.0],
##                 [-0.0 -0.0 -0.0],
##                 [-0.0 -0.0 -0.0],
##                 [-0.0 -0.0 -0.0]]

        assert num.allclose(num.min(w), -517.877771593)
        assert num.allclose(num.max(w), 0.0)
        assert num.allclose(slide.a3D, 518.38797486)

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import set_quantity [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_quantity [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_quantity [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_quantity [as 别名]
def setup_domain(simulation):
    
    args = simulation.args
    verbose = args.verbose
    alg = args.alg
    
    N = args.N
    S = args.S
    E = args.E
    W = args.W
    
    from catchment_info import create_catchment_list
    from catchment_info import create_manning_list
    
    CatchmentList = create_catchment_list(simulation)
    ManningList = create_manning_list(simulation)
    
    #------------------------------------------------------------------------------
    # CREATING MESH
    #------------------------------------------------------------------------------
    
    bounding_polygon = [[W, S], [E, S], [E, N], [W, N]]
    #interior_regions = read_polygon_dir(CatchmentDictionary, join('Model', 'Bdy'))
    interior_regions = read_polygon_list(CatchmentList)

    # FIXME: Have these in a shapefile / other file and read them in    
    breaklines=[[[306612.336559443,6193708.75358065],
                 [306604.441364239,6193693.17994946]],
                [[306977.886673843,6193753.44134088],
                 [306978.027867398,6193710.94208076]],
                [[306956.001672788,6193750.89985688],
                 [306956.707640564,6193706.14149989]],
                [[306627.303076293,6193697.45809624],
                 [306620.525785644,6193683.62112783]],
                [[307236.83565407,6193741.01630802],
                 [307231.682089306,6193721.03741996]],
                [[307224.975395434,6193742.71063068],
                 [307220.880782334,6193723.36711362]],
                [[307624.764946969,6193615.98941489],
                 [307617.98765632,6193601.44647871]],
                [[307613.328268998,6193623.19028621],
                 [307607.751123568,6193610.97704368]]]

    # Make the mesh
    create_mesh_from_regions(bounding_polygon, 
        boundary_tags={'south': [0], 'east': [1], 'north': [2], 'west': [3]},
        maximum_triangle_area=args.maximum_triangle_area,
        interior_regions=interior_regions,
        filename=args.meshname,
        breaklines=breaklines,
        use_cache=False,
        verbose=True)
    
    #------------------------------------------------------------------------------
    # SETUP COMPUTATIONAL DOMAIN
    #------------------------------------------------------------------------------
    
    domain = Domain(args.meshname, use_cache=False, verbose=True)

    domain.set_flow_algorithm(alg)

    if(not domain.get_using_discontinuous_elevation()):
        raise Exception, 'This model run relies on a discontinuous elevation solver (because of how topography is set up)'

    domain.set_datadir(args.model_output_dir)
    domain.set_name(args.outname)
        
    print domain.statistics()
    
    #------------------------------------------------------------------------------
    # APPLY MANNING'S ROUGHNESSES
    #------------------------------------------------------------------------------
    
    if verbose: print 'Calculating complicated polygon friction function'
    friction_list = read_polygon_list(ManningList)
    domain.set_quantity('friction', Polygon_function(friction_list, default=args.base_friction, geo_reference=domain.geo_reference))
    
    # Set a Initial Water Level over the Domain
    domain.set_quantity('stage', 0)
   
    # Decompress the zip file to make a csv for reading 
    zipfile.ZipFile('DEM_bridges/towradgi_cleaner.zip').extract('towradgi.csv',path='DEM_bridges/')

    if verbose: print 'Setting up elevation interpolation function'
    from anuga.utilities.quantity_setting_functions import make_nearestNeighbour_quantity_function
    elev_xyz=numpy.genfromtxt(fname=args.basename+'.csv',delimiter=',')

    # Use nearest-neighbour interpolation of elevation 
    elev_fun_wrapper=make_nearestNeighbour_quantity_function(elev_xyz,domain)
    if verbose: print 'Applying elevation interpolation function'    
    domain.set_quantity('elevation', elev_fun_wrapper, location='centroids')

    os.remove('DEM_bridges/towradgi.csv') # Clean up csv file
    
    return domain

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import set_quantity [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_quantity [as 别名]
domain.set_name(output_file)                
domain.set_datadir(output_dir) 

#------------------------------------------------------------------------------
# Setup Algorithm, either using command line arguments
# or override manually yourself
#------------------------------------------------------------------------------
from anuga.utilities.argparsing import parse_standard_args
alg = parse_standard_args()
domain.set_flow_algorithm(alg)


#------------------------------------------------------------------------------
# Setup initial conditions
#------------------------------------------------------------------------------
domain.set_quantity('elevation',0.0)
domain.set_quantity('friction', 0.0)

h0 = 5.0
h1 = 1.0

def height(x,y):
    z = zeros(len(x), float)
    for i in range(len(x)):
        if x[i]<=0.0:
            z[i] = h0
        else:
            z[i] = h1
    return z
domain.set_quantity('stage', height)

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import set_quantity [as 别名]
    
    domain = Domain(points, vertices, boundary) 
    
    domain.set_name(output_file)                
    domain.set_datadir(output_dir)  
    
    #------------------------------------------------------------------------------
    # Setup Algorithm, either using command line arguments
    # or override manually yourself
    #------------------------------------------------------------------------------
    domain.set_flow_algorithm(alg)
    
    #------------------------------------------------------------------------------
    # Setup initial conditions
    #------------------------------------------------------------------------------
    domain.set_quantity('elevation',-100.0)
    domain.set_quantity('friction', 0.00)
    domain.set_quantity('stage', 0.0)            

else: 

    domain = None
    
#-----------------------------------------------------------------------------------
# Parallel Domain 
#-----------------------------------------------------------------------------------
domain = distribute(domain)
   
#-----------------------------------------------------------------------------
# Setup boundary conditions
#------------------------------------------------------------------------------

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import set_quantity [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_quantity [as 别名]
evolved_quantities = ['stage', 'xmomentum', 'ymomentum', 'elevation']
                                               
domain = Domain(points, vertices, boundary, evolved_quantities=evolved_quantities)
domain.set_flow_algorithm('DE0')
domain.set_name('veg') # Output name
# domain.set_store_vertices_uniquely(True)

# print domain.statistics()

domain.set_quantities_to_be_stored({'elevation': 2,
                                    'stage': 2, 
                                    'xmomentum': 2,
                                    'ymomentum': 2})


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

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import set_quantity [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_quantity [as 别名]
    def test_erosion_operator_simple_1_5(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)
        domain.set_flow_algorithm('1_5')
        
        #Flat surface with 1m of water
        domain.set_quantity('elevation', 0.5)
        domain.set_quantity('stage', 1.0)
        domain.set_quantity('friction', 0)
        domain.set_quantity('xmomentum',2.0)
        domain.set_quantity('ymomentum',3.0)

        Stage = domain.quantities['stage'].centroid_values
        Elevation = domain.quantities['elevation'].centroid_values

        Height = Stage - Elevation

        sum1 = num.sum(Height)

        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]


        operator = Erosion_operator(domain, indices=indices, logging=True)

        # Apply Operator
        domain.timestep = 2.0
        operator()

        elev_ex  = [ 0.05555556,  0.11111111,  0.27777778,  0.05555556]
        stage_ex = [ 0.55555556,  0.61111111,  0.77777778,  0.55555556]

        Stage = domain.quantities['stage'].centroid_values
        Elevation = domain.quantities['elevation'].centroid_values

        Height = Stage - Elevation

        sum2 = num.sum(Height)
        
        #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 sum1 == sum2
        assert num.allclose(domain.quantities['stage'].centroid_values, stage_ex)
        assert num.allclose(domain.quantities['xmomentum'].centroid_values, 2.0)
        assert num.allclose(domain.quantities['ymomentum'].centroid_values, 3.0)

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import set_quantity [as 别名]
    domain.set_datadir(output_dir) 
    domain.set_flow_algorithm(alg)

    
    #------------------------------------------------------------------------------
    # Setup initial conditions
    #------------------------------------------------------------------------------
    def elevation(x,y):
        z_b = zeros(len(x))
        for i in range(len(x)):
            if (8.0 <= x[i] <= 12.0):
                z_b[i] = 0.2 - 0.05*(x[i]-10.0)**2
            else:
                z_b[i] = 0.0
        return z_b
    domain.set_quantity('elevation',elevation)
    domain.set_quantity('friction', 0.0)
    domain.set_quantity('xmomentum', 0.0)
    
    
    def stage(x,y):
        return 0.66*ones(len(x))
    domain.set_quantity('stage', stage)
else:
    domain = None
    
#-----------------------------------------------------------------------------
# create Parallel Domain
#------------------------------------------------------------------------------
domain = distribute(domain)