Python anuga.Domain类代码示例

    def test_create_operator(self):
        points = num.array([[0.0,0.0],[1.0,0.0],[0.0,1.0]])
        
        elements = num.array([[0,1,2]])
        boundary_map = {}
        boundary_map[(0,0)] = 'edge0'
        boundary_map[(0,1)] = 'edge1'
        boundary_map[(0,2)] = 'edge2'

        domain = Domain(points, elements, boundary_map)

        operator = Operator(domain)

        message = operator.statistics()
        assert message == 'You need to implement operator statistics for your operator'

        message = operator.timestepping_statistics()
        assert message == 'You need to implement timestepping statistics for your operator'

        domain.timestep = 3.0

        assert operator.get_timestep() == domain.get_timestep()
        
        try:
            operator()
        except:
            pass
        else:
            raise Exception('should have raised an exception')

    def __init__(self,
                 coordinates,
                 vertices,
                 boundary = None,
                 full_send_dict = None,
                 ghost_recv_dict = None,
                 velocity = None):

        Domain.__init__(self,
                        coordinates,
                        vertices,
                        boundary,
                        velocity = velocity,
                        full_send_dict=full_send_dict,
                        ghost_recv_dict=ghost_recv_dict,
                        processor=pypar.rank(),
                        numproc=pypar.size()
                        )

        N = self.number_of_elements


        self.communication_time = 0.0
        self.communication_reduce_time = 0.0


        print 'processor',self.processor
        print 'numproc',self.numproc

    def update_timestep(self, yieldstep, finaltime):
        """Calculate local timestep
        """

        generic_comms.communicate_flux_timestep(self, yieldstep, finaltime)

        Domain.update_timestep(self, yieldstep, finaltime)

    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)

    def set_name(self, name):
        """Assign name based on processor number 
        """

        if name.endswith('.sww'):
            name = name[:-4]

        self.global_name = name

        # Call parents method with processor number attached.
        Domain.set_name(self, name + '_P%d_%d' %(self.numproc, self.processor))

    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())

    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)

    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)

    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)

    def update_timestep(self, yieldstep, finaltime):

        #LINDA:
        # moved the calculation so that it is done after timestep
        # has been broadcast
        
#        # Calculate local timestep
#        Domain.update_timestep(self, yieldstep, finaltime)

        import time
        t0 = time.time()

        # For some reason it looks like pypar only reduces numeric arrays
        # hence we need to create some dummy arrays for communication
        ltimestep = num.ones( 1, num.float )
        ltimestep[0] = self.flux_timestep
        gtimestep = num.zeros( 1, num.float ) # Buffer for results

        #ltimestep = self.flux_timeste

        #print self.processor, ltimestep, gtimestep
        
        gtimestep = pypar.reduce(ltimestep, pypar.MIN, 0, buffer=gtimestep)

        #print self.processor, ltimestep, gtimestep
        
        pypar.broadcast(gtimestep,0)

        #print self.processor, ltimestep, gtimestep

        self.flux_timestep = gtimestep[0]
        
        self.communication_reduce_time += time.time()-t0

        # LINDA:
        # Now update time stats
        
        # Calculate local timestep
        Domain.update_timestep(self, yieldstep, finaltime)

    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 Reflective_boundary
from anuga import Dirichlet_boundary
from anuga import Time_boundary
import os


#------------------------------------------------------------------------------
# Setup computational domain
#------------------------------------------------------------------------------
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_name() # Output name based on script name. You can add timestamp=True
print domain.statistics()


#------------------------------------------------------------------------------
# Setup initial conditions
#------------------------------------------------------------------------------
def topography(x,y):
    """Complex topography defined by a function of vectors x and y."""

    z = -x/100

    # Step
    id = (2 < x) & (x < 4)
    z[id] += 0.4 - 0.05*y[id]

from anuga import rectangular_cross
from anuga import Domain
from anuga import Reflective_boundary
from anuga import Dirichlet_boundary
from anuga import Time_boundary

#------------------------------------------------------------------------------
# Setup computational domain
#------------------------------------------------------------------------------
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_name()
print domain.statistics()
domain.set_quantities_to_be_stored({'elevation': 2,
                                    'stage': 2,
                                    'xmomentum': 2,
                                    'ymomentum': 2})

#------------------------------------------------------------------------------
# Setup initial conditions
#------------------------------------------------------------------------------
def topography_dam(x,y):
    """Complex topography defined by a function of vectors x and y."""

    z = -x/100

interactive_visualisation = False

if myid == 0:
    #------------------------------------------------------------------------------
    # Setup sequential domain
    #------------------------------------------------------------------------------
    dx = 500.
    dy = dx
    L = 100000.
    W = 10*dx
    
    # structured mesh
    points, vertices, boundary = anuga.rectangular_cross(int(L/dx), int(W/dy), L, W, (0.0, -W/2))
    
    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)

    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)