Python Domain.set_quantities_to_be_stored方法代码示例

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import set_quantities_to_be_stored [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_quantities_to_be_stored [as 别名]
# Setup computational domain
#------------------------------------------------------------------------------
print ' Set up Domain first...'
length = 24.
width = 5.
dx = dy = 0.2 #.1           # Resolution: Length of subdivisions on both axes

points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
                                               len1=length, len2=width)
domain = Domain(points, vertices, boundary)
domain.set_flow_algorithm('DE1')
domain.set_name('flat_fill_slice_erosion') # Output name
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', expression='elevation')   # Dry initial condition

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

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import set_quantities_to_be_stored [as 别名]
from anuga import myid, finalize, distribute


args = anuga.get_args()
alg = args.alg
verbose = args.verbose

if myid == 0:
    # ---------
    # Setup computational domain
    # ---------
    points, vertices, boundary = anuga.rectangular_cross(100, 3, len1=1.0, len2=0.03)
    domain = Domain(points, vertices, boundary)  # Create Domain
    domain.set_name("runup")  # 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_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.

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import set_quantities_to_be_stored [as 别名]
evolved_quantities = ["stage", "xmomentum", "ymomentum", "elevation", "concentration"]

domain = Domain(points, vertices, boundary, evolved_quantities=evolved_quantities)
domain.set_flow_algorithm("DE0")
domain.set_name("test_equations")  # 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,
        "concentration": 2,
    }
)

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

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

#.........这里部分代码省略.........
            urs2sts(base_name,
                    basename_out=sts_file,
                    ordering_filename=order_file,
                    mean_stage=tide,
                    verbose=verbose)
            self.delete_mux(files)

            assert(os.access(sts_file+'.sts', os.F_OK))

            os.remove(order_file)

        barrier()
        boundary_polygon = create_sts_boundary(sts_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])

        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.
        if myid==0:
            create_mesh_from_regions(boundary_polygon,
                                     boundary_tags=boundary_tags,
                                     maximum_triangle_area=extent_res,
                                     filename=meshname,
                                     interior_regions=interior_regions,
                                     verbose=verbose)

        barrier()
        
        domain_fbound = Domain(meshname)
        domain_fbound.set_quantities_to_be_stored(None)
        domain_fbound.set_quantity('stage', tide)
        if verbose: print "Creating file boundary condition"
        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})

        temp_fbound=num.zeros(int(finaltime/yieldstep)+1,num.float)
        if verbose: print "Evolving domain with file boundary condition"
        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]
            if verbose: domain_fbound.write_time()
            
        
        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, f=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,f=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})
        temp_drchlt=num.zeros(int(finaltime/yieldstep)+1,num.float)
        
        for i, t in enumerate(domain_drchlt.evolve(yieldstep=yieldstep,
                                                   finaltime=finaltime, 
                                                   skip_initial_step = False)):
            temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2]
            #domain_drchlt.write_time()
        
        #print domain_fbound.quantities['stage'].vertex_values
        #print domain_drchlt.quantities['stage'].vertex_values
                    
        assert num.allclose(temp_fbound,temp_drchlt),temp_fbound-temp_drchlt

        
        assert num.allclose(domain_fbound.quantities['stage'].vertex_values,
                            domain_drchlt.quantities['stage'].vertex_values)
                        
        assert num.allclose(domain_fbound.quantities['xmomentum'].vertex_values,
                            domain_drchlt.quantities['xmomentum'].vertex_values)                        
                        
        assert num.allclose(domain_fbound.quantities['ymomentum'].vertex_values,
                            domain_drchlt.quantities['ymomentum'].vertex_values)
        
        if not sys.platform == 'win32':
            if myid==0: os.remove(sts_file+'.sts')
        
        if myid==0: os.remove(meshname)

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

#.........这里部分代码省略.........
        # 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])


        assert num.allclose(bounding_polygon_utm,boundary_polygon)

        extent_res=10000
        meshname = 'urs_test_mesh' + '.tsh'
        interior_regions=None
        boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]}
        
        #------------------------------------------------------------
        # Create mesh on the master processor and store in file. This file
        # is read in by each slave processor when needed
        #------------------------------------------------------------
        if myid==0:
            create_mesh_from_regions(boundary_polygon,
                                     boundary_tags=boundary_tags,
                                     maximum_triangle_area=extent_res,
                                     filename=meshname,
                                     interior_regions=interior_regions,
                                     verbose=verbose)
        

            # barrier()
            domain_fbound = Domain(meshname)
            domain_fbound.set_quantities_to_be_stored(None)
            domain_fbound.set_quantity('stage', tide)
            # print domain_fbound.mesh.get_boundary_polygon()
        else:
            domain_fbound=None

        barrier()
        if ( verbose and myid == 0 ): 
            print 'DISTRIBUTING PARALLEL DOMAIN'
        domain_fbound = distribute(domain_fbound)

        #--------------------------------------------------------------------
        # Find which sub_domain in which the interpolation points are located 
        #
        # Sometimes the interpolation points sit exactly
        # between two centroids, so in the parallel run we
        # reset the interpolation points to the centroids
        # found in the sequential run
        #--------------------------------------------------------------------
        interpolation_points = [[279000,664000], [280250,664130], 
                                    [279280,665400], [280500,665000]]

        interpolation_points=num.array(interpolation_points)

        #if myid==0:
        #    import pylab as P
        #    boundary_polygon=num.array(boundary_polygon)
        #    P.plot(boundary_polygon[:,0],boundary_polygon[:,1])
        #    P.plot(interpolation_points[:,0],interpolation_points[:,1],'ko')
        #    P.show()

        fbound_gauge_values = []
        fbound_proc_tri_ids = []

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

domain.set_flow_algorithm('1_75')
domain.set_name('run_simple_sed_transport_veg') # Output name
domain.set_store_vertices_uniquely(True)
domain.set_quantity('elevation', topography)           # elevation is a function
domain.set_quantity('stage', expression='elevation')   # Dry initial condition

print domain.statistics()

"""
Store process-specific quantities with same functions
""" 
domain.set_quantities_to_be_stored({'elevation': 2,
                                    'stage': 2,
                                    'xmomentum': 2,
                                    'ymomentum': 2,
                                    'concentration': 2,
                                    'vegetation': 1})
                                    
#------------------------------------------------------------------------------
# Setup boundary conditions
#------------------------------------------------------------------------------
max_elev = domain.quantities['elevation'].vertex_values.max()
min_elev = domain.quantities['elevation'].vertex_values.min()
"""
Use operator-specific Reflective and Dirichlet boundaries. Use original
Dirichlet boundary for outlet.

""" 
Bi = Dirichlet_boundary_Sed([max_elev + 0.5, 0, 0, 0.2])  # Inflow, 20% sed
Br = Reflective_boundary_Sed(domain)           # Solid reflective wall