Python Domain.timestep方法代码示例

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import timestep [as 别名]
    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')

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import timestep [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 timestep [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 timestep [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 timestep [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 timestep [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 timestep [as 别名]
    def test_rate_operator_functions_empty_indices(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 = []
        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

        domain.tri_full_flag[0] = 0
        operator = Rate_operator(domain, rate=main_spatial_rate, factor=factor, \
                      indices=indices, default_rate = default_rate)


        # Apply Operator
        domain.timestep = 2.0
        operator()

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

        #print Q_ex, Q
        #print indices
        #print "d"
        #print d
        stage_ex = num.array([ 1.0,  1.0,   1.0,  1.0])
        stage_ex[indices] = 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)
        assert num.allclose(domain.fractional_step_volume_integral, ((d-1.)*domain.areas[indices]).sum())

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

#.........这里部分代码省略.........

        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, \
                      indices=indices, default_rate = default_rate)


        # Apply Operator
        domain.set_time(360.0)
        domain.timestep = 1.0

        operator()



        d = domain.get_time()**2 * factor + 1.0
        stage_ex0 = [ d,  d,   1.,  d]

#        print d, domain.get_time(), F(360.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_ex0)
        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_time(-10.0)
        domain.timestep = 1.0

        try:
            operator()
        except:
            pass
        else:
            raise Exception('Should have raised an exception, time too early')


        domain.set_time(1300.0)
        domain.timestep = 1.0

        operator()

        d = default_rate*factor + d
        stage_ex1 = [ d,  d,   1.,  d]

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

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