Python Domain.get_quantity方法代码示例

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import get_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 get_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 get_quantity [as 别名]
#------------------------------------------------------------------------------
# Evolve system through time
#------------------------------------------------------------------------------

#from anuga.operators.set_elevation_operators import Circular_set_elevation_operator

#op1 = Circular_set_elevation_operator(domain, elevation=pole, radius=0.5, center = (12.0,3.0))


dam_break = False

for t in domain.evolve(yieldstep=0.1, finaltime=40.0):
    domain.print_timestepping_statistics()
    domain.print_operator_timestepping_statistics()

    if t >= 10 and not dam_break:
        print 'changing elevation'

        stage_c = domain.get_quantity('stage').centroid_values
        elev_c =  domain.get_quantity('elevation').centroid_values
        height_c = stage_c - elev_c
        domain.set_quantity('elevation', topography_dam_break)
        stage_c[:] = elev_c + height_c

        dam_break = True

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import get_quantity [as 别名]
#------------------------------------------------------------------------------
polygon1 = [ [10.0, 0.0], [11.0, 0.0], [11.0, 5.0], [10.0, 5.0] ]
polygon2 = [ [12.0, 2.0], [13.0, 2.0], [13.0, 3.0], [12.0, 3.0] ]

from anuga.operators.rate_operators import Rate_operator

op1 = Rate_operator(domain, rate=lambda t: 10.0 if (t>=0.0) else 0.0, polygon=polygon2)
op2 = Rate_operator(domain, rate=lambda t: 10.0 if (t>=0.0) else 0.0, radius=0.5, center=(10.0, 3.0))


domain.set_starttime(-0.1)
for t in domain.evolve(yieldstep=0.01, finaltime=0.0):
    domain.print_timestepping_statistics()
    domain.print_operator_timestepping_statistics()

    stage = domain.get_quantity('stage')
    elev  = domain.get_quantity('elevation')
    height = stage - elev

    print 'integral = ', height.get_integral()


for t in domain.evolve(yieldstep=0.1, duration=5.0):

    domain.print_timestepping_statistics()
    domain.print_operator_timestepping_statistics()

    stage = domain.get_quantity('stage')
    elev  = domain.get_quantity('elevation')
    height = stage - elev

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import get_quantity [as 别名]
    domain.set_datadir(".")  # Use current folder
    domain.set_quantities_to_be_stored({"stage": 2, "xmomentum": 2, "ymomentum": 2, "elevation": 1})
    domain.set_flow_algorithm(alg)

    # ------------------
    # Define topography
    # ------------------
    def topography(x, y):
        return -x / 2  # Linear bed slope

    def stagefun(x, y):
        return -0.45  # Stage

    domain.set_quantity("elevation", topography)  # Use function for elevation
    domain.get_quantity(
        "elevation"
    ).smooth_vertex_values()  # Steve's fix -- without this, substantial artificial velcities are generated everywhere in the domain. With this fix, there are artificial velocities near the coast, but not elsewhere.
    domain.set_quantity("friction", 0.0)  # Constant friction
    domain.set_quantity("stage", stagefun)  # Constant negative initial stage
else:
    domain = None

# --------------------------
# create Parallel Domain
# --------------------------
domain = distribute(domain)

# Setup boundary conditions
# --------------------------
Br = anuga.Reflective_boundary(domain)  # Solid reflective wall
Bt = anuga.Transmissive_boundary(domain)  # Continue all values of boundary -- not used in this example

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

#.........这里部分代码省略.........
        fbound_gauge_values = []
        fbound_proc_tri_ids = []
        for i, point in enumerate(interpolation_points):
            fbound_gauge_values.append([]) # Empty list for timeseries

            try:
                k = domain_fbound.get_triangle_containing_point(point)
                if domain_fbound.tri_full_flag[k] == 1:
                    fbound_proc_tri_ids.append(k)
                else:
                    fbound_proc_tri_ids.append(-1)            
            except:
                fbound_proc_tri_ids.append(-2)


        if verbose: print 'P%d has points = %s' %(myid, fbound_proc_tri_ids)

        #------------------------------------------------------------
        # Set boundary conditions
        #------------------------------------------------------------
        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})

        #------------------------------------------------------------
        # Evolve the domain on each processor
        #------------------------------------------------------------  
        for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep,
                                                   finaltime=finaltime, 
                                                   skip_initial_step = False)):

            stage = domain_fbound.get_quantity('stage')
            for i in range(4):
                if fbound_proc_tri_ids[i] > -1:
                    fbound_gauge_values[i].append(stage.centroid_values[fbound_proc_tri_ids[i]])
        
        #------------------------------------------------------------
        # Create domain to be run sequntially on each processor
        #------------------------------------------------------------
        domain_drchlt = Domain(meshname)
        domain_drchlt.set_quantities_to_be_stored(None)
        domain_drchlt.set_starttime(time_step)
        domain_drchlt.set_quantity('stage', tide)
        Br = Reflective_boundary(domain_drchlt)
        #Bd = Dirichlet_boundary([2.0+tide,220+10*tide,-220-10*tide])
        Bd = Time_boundary(domain=domain_drchlt, function=lambda t: [2.0+t/finaltime+tide,220.+10.*tide+10.*t/finaltime,-220.-10.*tide-10.*t/finaltime])
        #Bd = Time_boundary(domain=domain_drchlt,function=lambda t: [2.0+num.sin(t)+tide,10.*(2+20.+num.sin(t)+tide),-10.*(2+20.+num.sin(t)+tide)])
        domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br})
       
        drchlt_gauge_values = []
        drchlt_proc_tri_ids = []
        for i, point in enumerate(interpolation_points):
            drchlt_gauge_values.append([]) # Empty list for timeseries

            try:
                k = domain_drchlt.get_triangle_containing_point(point)
                if domain_drchlt.tri_full_flag[k] == 1:
                    drchlt_proc_tri_ids.append(k)
                else:
                    drchlt_proc_tri_ids.append(-1)            
            except:
                drchlt_proc_tri_ids.append(-2)


        if verbose: print 'P%d has points = %s' %(myid, drchlt_proc_tri_ids)

        #------------------------------------------------------------
        # Evolve entire domain on each processor
        #------------------------------------------------------------
        for i, t in enumerate(domain_drchlt.evolve(yieldstep=yieldstep,
                                                   finaltime=finaltime, 
                                                   skip_initial_step = False)):

            stage = domain_drchlt.get_quantity('stage')
            for i in range(4):
                drchlt_gauge_values[i].append(stage.centroid_values[drchlt_proc_tri_ids[i]])

        #------------------------------------------------------------
        # Compare sequential values with parallel values
        #------------------------------------------------------------
        barrier()
        success = True
        for i in range(4):
            if fbound_proc_tri_ids[i] > -1:
                fbound_gauge_values[i]=num.array(fbound_gauge_values[i])
                drchlt_gauge_values[i]=num.array(drchlt_gauge_values[i])
                #print i,fbound_gauge_values[i][4]
                #print i,drchlt_gauge_values[i][4]
                success = success and num.allclose(fbound_gauge_values[i], drchlt_gauge_values[i])
                assert success#, (fbound_gauge_values[i]-drchlt_gauge_values[i])

        #assert_(success)       

        if not sys.platform == 'win32':
            if myid==0: os.remove(sts_file+'.sts')
        
        if myid==0: os.remove(meshname)