Changeset 229


Ignore:
Timestamp:
Aug 26, 2004, 10:14:19 PM (21 years ago)
Author:
ole
Message:

Finally quashed bug where ymomentum wasn't updated in
second order limiter

Location:
inundation/ga/storm_surge/pyvolution
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • inundation/ga/storm_surge/pyvolution/domain.py

    r218 r229  
    290290            #Update boundary values
    291291            self.update_boundary()           
     292           
     293            #print
     294            #for name in self.conserved_quantities:
     295            #    Q = self.quantities[name]
     296            #    print 'Vertices (%s):' %name, Q.vertex_values[:4]                         
     297           
     298            #print
     299            #for name in self.conserved_quantities:
     300            #    Q = self.quantities[name]
     301            #    print 'Edges (%s):' %name, Q.edge_values[:4]               
    292302
    293303            #Compute all fluxes and timestep suitable for all volumes
    294304            self.compute_fluxes()
     305            ##print
     306            ##for name in self.conserved_quantities:
     307            ##    Q = self.quantities[name]
     308            ##    print 'EU:', Q.explicit_update[:4]       
    295309
    296310            #Update timestep to fit yieldstep and finaltime
    … …  
    300314            self.update_conserved_quantities()           
    301315
     316            #print
     317            #print 'Centroids'
     318            #print self.quantities['level'].centroid_values[1:4],\
     319            #      self.quantities['level'].centroid_values[13]
     320            #print self.quantities['xmomentum'].centroid_values[1:4],\
     321            #      self.quantities['xmomentum'].centroid_values[13]         
     322            #print self.quantities['ymomentum'].centroid_values[1:4],\
     323            #      self.quantities['ymomentum'].centroid_values[13]
     324                                                           
     325           
    302326            #Update vertex and edge values
    303327            self.distribute_to_vertices_and_edges()
     328           
     329            #print
     330            #for name in self.conserved_quantities:
     331            #    Q = self.quantities[name]
     332            #    print 'Vertices (%s):' %name, Q.vertex_values[1:4], Q.vertex_values[13]
     333                                           
     334           
     335           
    304336               
    305337            #Update time   
    … …  
    310342                self.number_of_first_order_steps += 1       
    311343
     344               
     345            #print 'Time=', self.time, self.timestep   
     346            #print self.quantities['level'].centroid_values[:4]
     347            #print self.quantities['xmomentum'].centroid_values[:4]
     348            #print self.quantities['ymomentum'].centroid_values[:4]                             
     349            #print
     350               
    312351            #Yield results
    313352            if finaltime is not None and abs(self.time - finaltime) < epsilon:

  • inundation/ga/storm_surge/pyvolution/quantity.py

    r221 r229  
    7878            v2 = self.vertex_values[i, 2]
    7979           
    80             self.centroid_values[i] = (v0 + v1 + v2)/3.0
    81            
     80            self.centroid_values[i] = (v0 + v1 + v2)/3
     81           
     82            #FIXME: This is a duplicate from 'interpolate_from_vertices_to_edges'
    8283            self.edge_values[i, 0] = 0.5*(v1 + v2)
    8384            self.edge_values[i, 1] = 0.5*(v0 + v2)           

  • inundation/ga/storm_surge/pyvolution/shallow_water.py

    r221 r229  
    273273def distribute_to_vertices_and_edges(domain):
    274274
    275     #print 'Calling distrib within SW'
     275    ##print 'Calling distrib within SW'
    276276    if domain.order == 1:
     277        #FIXME: This should be cleaned up, but we try to follow
     278        #pyvolution 2 strictly for now
     279        protect_against_negative_heights_centroid(domain)       
    277280        protect_against_infinitesimal_heights_centroid(domain)       
    278281        extrapolate_first_order(domain)
    279282    elif domain.order == 2:
    280         protect_against_negative_heights_centroid(domain)       
     283        #protect_against_negative_heights_centroid(domain)       
    281284        extrapolate_second_order(domain)       
    282285    else:
    … …  
    323326    """
    324327
    325     #FIXME: USed only in second order
     328   
    326329    #Water levels at centroids
    327330    wc = domain.quantities['level'].centroid_values
    … …  
    428431        #Update water levels
    429432        for i in range(3):
    430             #wv[k,i] = zv[k,i] + hc[k] + alpha*(zc[k]-zv[k,i])
    431             wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k])
    432 
    433 
     433            #FIXME: Use the original first-order one first, then switch
     434            wv[k,i] = zv[k,i] + hc[k] + alpha*(zc[k]-zv[k,i])
     435            #wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k])
     436
     437        #FIXME: What about alpha weighting of momentum??
     438           
    434439       
    435440def extrapolate_second_order(domain):
    … …  
    462467        Q = domain.quantities[name]
    463468        Q.extrapolate_second_order()
     469       
     470    #print 'y1', Q.vertex_values[1,:]   #OK
     471   
     472    #FIXME - like pyvolution 2 .......................
     473    protect_against_negative_heights_centroid(domain)   
     474
     475    #print 'y1', Q.vertex_values[1,:]   #OK   
     476   
     477    for name in domain.conserved_quantities:
     478        Q = domain.quantities[name]     
    464479        Q.limit()
    465480
    466481 
     482    #print 'y1', Q.vertex_values[1,:]   #OK   
     483     
    467484    #Water levels at centroids
    468485    wc = domain.quantities['level'].centroid_values
    … …  
    507524        #first order scheme and alpha==1 means using the limited
    508525        #second order scheme for the stage w
    509         #If hmin < 0 then alpha=0 reverrting to first order.
     526        #If hmin < 0 then alpha=0 reverting to first order.
    510527     
    511528        if z_range > 0.0:
    … …  
    513530        else:
    514531            alpha = 1.0
    515 
     532   
     533        ##if k==1: print 'alpha', alpha #OK
     534       
    516535        #Update water levels
    517536        #  (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) =
    518         #   zvi + hc + alpha*(hvi - hc)         
    519         for i in range(3):
    520             wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k])
    521            
    522 
    523         #Momentums at centroids
    524         xmomc = domain.quantities['xmomentum'].centroid_values
    525         ymomc = domain.quantities['ymomentum'].centroid_values       
    526 
    527         #Momentums at vertices
    528         xmomv = domain.quantities['xmomentum'].vertex_values
    529         ymomv = domain.quantities['ymomentum'].vertex_values               
    530 
    531         # Update momentum as a linear combination of
    532         # xmomc and ymomc (shallow) and momentum
    533         # from extrapolator xmomv and ymomv (deep).
    534        
    535         xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:];           
    536            
    537 
     537        #   zvi + hc + alpha*(hvi - hc)
     538        if alpha < 1:         
     539            for i in range(3):
     540                wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k])
     541           
     542
     543            #Momentums at centroids
     544            xmomc = domain.quantities['xmomentum'].centroid_values
     545            ymomc = domain.quantities['ymomentum'].centroid_values       
     546
     547            #Momentums at vertices
     548            xmomv = domain.quantities['xmomentum'].vertex_values
     549            ymomv = domain.quantities['ymomentum'].vertex_values         
     550
     551            # Update momentum as a linear combination of
     552            # xmomc and ymomc (shallow) and momentum
     553            # from extrapolator xmomv and ymomv (deep).
     554       
     555            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:];           
     556            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:];                     
     557           
     558    #print 'y1', Q.vertex_values[1,:]   #OK   
     559           
    538560    #Finally, protect against infinitesimal heights and high speeds
    539561    #Water depths at vertices

  • inundation/ga/storm_surge/pyvolution/test_shallow_water.py

    r221 r229  
    859859        assert allclose(L[1], [4.1, 16.1, 20.1])
    860860
    861 #FIXME: Make version with variable height!!!!!!!!
     861 
     862
     863    def test_second_order_distribute_real_data(self):
     864        #Using test data generated by pyvolution-2
     865        #Assuming no friction and flat bed (0.0)
     866
     867        a = [0.0, 0.0]
     868        b = [0.0, 1.0/5]
     869        c = [0.0, 2.0/5]
     870        d = [1.0/5, 0.0]
     871        e = [1.0/5, 1.0/5]
     872        f = [1.0/5, 2.0/5]
     873        g = [2.0/5, 2.0/5]
     874
     875        points = [a, b, c, d, e, f, g]
     876        #bae, efb, cbf, feg
     877        vertices = [ [1,0,4], [4,5,1], [2,1,5], [5,4,6]]
     878       
     879        domain = Domain(points, vertices)
     880
     881        def slope(x, y):
     882            return -x/3
     883
     884        domain.set_quantity('elevation', slope)
     885        domain.set_quantity('level',
     886                            [0.01298164, 0.00365611, 0.01440365, -0.0381856437096],
     887                            'centroids')
     888        domain.set_quantity('xmomentum',
     889                            [0.00670439, 0.01263789, 0.00647805, 0.0178180740668],
     890                            'centroids')
     891        domain.set_quantity('ymomentum',                           
     892                            [-7.23510980e-004, -6.30413883e-005, 6.30413883e-005, 0.000200907255866],
     893                            'centroids')                           
     894       
     895        E = domain.quantities['elevation'].vertex_values
     896        L = domain.quantities['level'].vertex_values
     897        X = domain.quantities['xmomentum'].vertex_values       
     898        Y = domain.quantities['ymomentum'].vertex_values               
     899
     900        #print E
     901        domain.order = 2       
     902        domain.distribute_to_vertices_and_edges()
     903
     904        #print L[1,:]
     905        #print X[1,:]
     906        #print Y[1,:]
     907
     908        assert allclose(L[1,:], [-0.00825735775384, -0.00801881482869, 0.0272445025825])
     909        assert allclose(X[1,:], [0.0143507718962, 0.0142502147066, 0.00931268339717])
     910        assert allclose(Y[1,:], [-0.000117062180693, 7.94434448109e-005, -0.000151505429018])
     911       
    862912       
    863913    def test_second_order_flat_bed_onestep(self):
    … …  
    12971347                         -0.14183568, -0.13623283, -0.14529803])
    12981348       
    1299     def test_bedslope_problem_second_order(self):
     1349    def test_bedslope_problem_second_order_one_step(self):
    13001350
    13011351        from mesh_factory import rectangular
    … …  
    13571407       
    13581408        #Evolution
    1359         for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
     1409        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05):
    13601410            #domain.write_time()           
    13611411            pass
    13621412
    13631413
    1364         #print domain.quantities['level'].centroid_values
    1365        
    1366         #Data from earlier version of pyvolution
     1414        #print domain.quantities['level'].centroid_values
     1415        assert allclose(domain.quantities['level'].centroid_values,
     1416                 [0.01290985, 0.02356019, 0.01619096, 0.02356019, 0.01619096,
     1417                  0.02356019, 0.01619096, 0.02356019, 0.01619096, 0.02356019,
     1418                  0.01619096, 0.0268413, -0.04411074, -0.0248011, -0.04186556,
     1419                  -0.0248011, -0.04186556, -0.0248011, -0.04186556, -0.0248011,
     1420                  -0.04186556, -0.0248011, -0.04186556, -0.02255593,
     1421                  -0.09966629, -0.08035666, -0.09742112, -0.08035666,
     1422                  -0.09742112, -0.08035666, -0.09742112, -0.08035666,
     1423                  -0.09742112, -0.08035666, -0.09742112, -0.07811149,
     1424                  -0.15522185, -0.13591222, -0.15297667, -0.13591222,
     1425                  -0.15297667, -0.13591222, -0.15297667, -0.13591222,
     1426                  -0.15297667, -0.13591222, -0.15297667, -0.13366704,
     1427                  -0.2107774, -0.19146777, -0.20853223, -0.19146777,
     1428                  -0.20853223, -0.19146777, -0.20853223, -0.19146777,
     1429                  -0.20853223, -0.19146777, -0.20853223, -0.1892226,
     1430                  -0.26120669, -0.24776246, -0.25865535, -0.24776246,
     1431                  -0.25865535, -0.24776246, -0.25865535, -0.24776246,
     1432                  -0.25865535, -0.24776246, -0.25865535, -0.24521113])
     1433
     1434       
     1435
     1436    def test_bedslope_problem_second_order_two_steps(self):
     1437
     1438        from mesh_factory import rectangular
     1439        from shallow_water import Domain, Reflective_boundary, Constant_height
     1440        from Numeric import array
     1441
     1442        #Create basic mesh
     1443        points, vertices, boundary = rectangular(6, 6)
     1444
     1445        #Create shallow water domain
     1446        domain = Domain(points, vertices, boundary)
     1447        domain.smooth = False
     1448        domain.default_order=2
     1449
     1450        #Bed-slope and friction at vertices (and interpolated elsewhere)
     1451        def x_slope(x, y):
     1452            return -x/3
     1453
     1454        domain.set_quantity('elevation', x_slope)
     1455
     1456        # Boundary conditions
     1457        Br = Reflective_boundary(domain)
     1458        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
     1459
     1460        #Initial condition
     1461        domain.set_quantity('level', Constant_height(x_slope, 0.05))
     1462        domain.check_integrity()
     1463
     1464        assert allclose(domain.quantities['level'].centroid_values,
     1465                        [0.01296296, 0.03148148, 0.01296296,
     1466                        0.03148148, 0.01296296, 0.03148148,
     1467                        0.01296296, 0.03148148, 0.01296296,
     1468                        0.03148148, 0.01296296, 0.03148148,
     1469                        -0.04259259, -0.02407407, -0.04259259,
     1470                        -0.02407407, -0.04259259, -0.02407407,
     1471                        -0.04259259, -0.02407407, -0.04259259,
     1472                        -0.02407407, -0.04259259, -0.02407407,
     1473                        -0.09814815, -0.07962963, -0.09814815,
     1474                        -0.07962963, -0.09814815, -0.07962963,
     1475                        -0.09814815, -0.07962963, -0.09814815,
     1476                        -0.07962963, -0.09814815, -0.07962963,
     1477                        -0.1537037 , -0.13518519, -0.1537037,
     1478                        -0.13518519, -0.1537037, -0.13518519,
     1479                        -0.1537037 , -0.13518519, -0.1537037,
     1480                        -0.13518519, -0.1537037, -0.13518519,
     1481                        -0.20925926, -0.19074074, -0.20925926,
     1482                        -0.19074074, -0.20925926, -0.19074074,
     1483                        -0.20925926, -0.19074074, -0.20925926,
     1484                        -0.19074074, -0.20925926, -0.19074074,
     1485                        -0.26481481, -0.2462963, -0.26481481,
     1486                        -0.2462963, -0.26481481, -0.2462963,
     1487                        -0.26481481, -0.2462963, -0.26481481,
     1488                        -0.2462963, -0.26481481, -0.2462963])
     1489
     1490
     1491        #print domain.quantities['level'].extrapolate_second_order()
     1492        #domain.distribute_to_vertices_and_edges()
     1493        #print domain.quantities['level'].vertex_values[:,0]       
     1494       
     1495        #Evolution
     1496        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
     1497            pass
     1498
     1499       
     1500        #Data from earlier version of pyvolution ft=0.1
    13671501        assert allclose(domain.min_timestep, 0.0376895634803)
    13681502        assert allclose(domain.max_timestep, 0.0415635655309)
    13691503
    13701504
    1371         #FIXME: Fails
    1372         #assert allclose(domain.quantities['level'].centroid_values,
    1373         #                [0.00855788, 0.01575204, 0.00994606, 0.01717072,
    1374         #                 0.01005985, 0.01716362, 0.01005985, 0.01716299,
    1375         #                 0.01007098, 0.01736248, 0.01216452, 0.02026776,
    1376         #                 -0.04462374, -0.02479045, -0.04199789, -0.0229465,
    1377         #                 -0.04184033, -0.02295693, -0.04184013, -0.02295675,
    1378         #                 -0.04184486, -0.0228168, -0.04028876, -0.02036486,
    1379         #                 -0.10029444, -0.08170809, -0.09772846, -0.08021704,
    1380         #                 -0.09760006, -0.08022143, -0.09759984, -0.08022124,
    1381         #                 -0.09760261, -0.08008893, -0.09603914, -0.07758209,
    1382         #                 -0.15584152, -0.13723138, -0.15327266, -0.13572906,
    1383         #                 -0.15314427, -0.13573349, -0.15314405, -0.13573331,
    1384         #                 -0.15314679, -0.13560104, -0.15158523, -0.13310701,
    1385         #                 -0.21208605, -0.19283913, -0.20955631, -0.19134189,
    1386         #                 -0.20942821, -0.19134598, -0.20942799, -0.1913458,
    1387         #                 -0.20943005, -0.19120952, -0.20781177, -0.18869401,
    1388         #                 -0.25384082, -0.2463294, -0.25047649, -0.24464654,
    1389         #                 -0.25031159, -0.24464253, -0.25031112, -0.24464253,
    1390         #                 -0.25031463, -0.24454764, -0.24885323, -0.24286438])
     1505        assert allclose(domain.quantities['level'].centroid_values,
     1506                        [0.00855788, 0.01575204, 0.00994606, 0.01717072,
     1507                         0.01005985, 0.01716362, 0.01005985, 0.01716299,
     1508                         0.01007098, 0.01736248, 0.01216452, 0.02026776,
     1509                         -0.04462374, -0.02479045, -0.04199789, -0.0229465,
     1510                         -0.04184033, -0.02295693, -0.04184013, -0.02295675,
     1511                         -0.04184486, -0.0228168, -0.04028876, -0.02036486,
     1512                         -0.10029444, -0.08170809, -0.09772846, -0.08021704,
     1513                         -0.09760006, -0.08022143, -0.09759984, -0.08022124,
     1514                         -0.09760261, -0.08008893, -0.09603914, -0.07758209,
     1515                         -0.15584152, -0.13723138, -0.15327266, -0.13572906,
     1516                         -0.15314427, -0.13573349, -0.15314405, -0.13573331,
     1517                         -0.15314679, -0.13560104, -0.15158523, -0.13310701,
     1518                         -0.21208605, -0.19283913, -0.20955631, -0.19134189,
     1519                         -0.20942821, -0.19134598, -0.20942799, -0.1913458,
     1520                         -0.20943005, -0.19120952, -0.20781177, -0.18869401,
     1521                         -0.25384082, -0.2463294, -0.25047649, -0.24464654,
     1522                         -0.25031159, -0.24464253, -0.25031112, -0.24464253,
     1523                         -0.25031463, -0.24454764, -0.24885323, -0.24286438])
     1524
     1525       
     1526       
     1527               
     1528    def test_bedslope_problem_second_order_two_yieldsteps(self):
     1529
     1530        from mesh_factory import rectangular
     1531        from shallow_water import Domain, Reflective_boundary, Constant_height
     1532        from Numeric import array
     1533
     1534        #Create basic mesh
     1535        points, vertices, boundary = rectangular(6, 6)
     1536
     1537        #Create shallow water domain
     1538        domain = Domain(points, vertices, boundary)
     1539        domain.smooth = False
     1540        domain.default_order=2
     1541
     1542        #Bed-slope and friction at vertices (and interpolated elsewhere)
     1543        def x_slope(x, y):
     1544            return -x/3
     1545
     1546        domain.set_quantity('elevation', x_slope)
     1547
     1548        # Boundary conditions
     1549        Br = Reflective_boundary(domain)
     1550        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
     1551
     1552        #Initial condition
     1553        domain.set_quantity('level', Constant_height(x_slope, 0.05))
     1554        domain.check_integrity()
     1555
     1556        assert allclose(domain.quantities['level'].centroid_values,
     1557                        [0.01296296, 0.03148148, 0.01296296,
     1558                        0.03148148, 0.01296296, 0.03148148,
     1559                        0.01296296, 0.03148148, 0.01296296,
     1560                        0.03148148, 0.01296296, 0.03148148,
     1561                        -0.04259259, -0.02407407, -0.04259259,
     1562                        -0.02407407, -0.04259259, -0.02407407,
     1563                        -0.04259259, -0.02407407, -0.04259259,
     1564                        -0.02407407, -0.04259259, -0.02407407,
     1565                        -0.09814815, -0.07962963, -0.09814815,
     1566                        -0.07962963, -0.09814815, -0.07962963,
     1567                        -0.09814815, -0.07962963, -0.09814815,
     1568                        -0.07962963, -0.09814815, -0.07962963,
     1569                        -0.1537037 , -0.13518519, -0.1537037,
     1570                        -0.13518519, -0.1537037, -0.13518519,
     1571                        -0.1537037 , -0.13518519, -0.1537037,
     1572                        -0.13518519, -0.1537037, -0.13518519,
     1573                        -0.20925926, -0.19074074, -0.20925926,
     1574                        -0.19074074, -0.20925926, -0.19074074,
     1575                        -0.20925926, -0.19074074, -0.20925926,
     1576                        -0.19074074, -0.20925926, -0.19074074,
     1577                        -0.26481481, -0.2462963, -0.26481481,
     1578                        -0.2462963, -0.26481481, -0.2462963,
     1579                        -0.26481481, -0.2462963, -0.26481481,
     1580                        -0.2462963, -0.26481481, -0.2462963])
     1581
     1582
     1583        #print domain.quantities['level'].extrapolate_second_order()
     1584        #domain.distribute_to_vertices_and_edges()
     1585        #print domain.quantities['level'].vertex_values[:,0]       
     1586       
     1587        #Evolution
     1588        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):   #0.05??
     1589            #domain.write_time()           
     1590            pass
     1591
     1592
     1593       
     1594        assert allclose(domain.quantities['level'].centroid_values,
     1595                 [0.00855788, 0.01575204, 0.00994606, 0.01717072, 0.01005985,
     1596                  0.01716362, 0.01005985, 0.01716299, 0.01007098, 0.01736248,
     1597                  0.01216452, 0.02026776, -0.04462374, -0.02479045, -0.04199789,
     1598                  -0.0229465, -0.04184033, -0.02295693, -0.04184013,
     1599                  -0.02295675, -0.04184486, -0.0228168, -0.04028876,
     1600                  -0.02036486, -0.10029444, -0.08170809, -0.09772846,
     1601                  -0.08021704, -0.09760006, -0.08022143, -0.09759984,
     1602                  -0.08022124, -0.09760261, -0.08008893, -0.09603914,
     1603                  -0.07758209, -0.15584152, -0.13723138, -0.15327266,
     1604                  -0.13572906, -0.15314427, -0.13573349, -0.15314405,
     1605                  -0.13573331, -0.15314679, -0.13560104, -0.15158523,
     1606                  -0.13310701, -0.21208605, -0.19283913, -0.20955631,
     1607                  -0.19134189, -0.20942821, -0.19134598, -0.20942799,
     1608                  -0.1913458, -0.20943005, -0.19120952, -0.20781177,
     1609                  -0.18869401, -0.25384082, -0.2463294, -0.25047649,
     1610                  -0.24464654, -0.25031159, -0.24464253, -0.25031112,
     1611                  -0.24464253, -0.25031463, -0.24454764, -0.24885323,
     1612                  -0.24286438])
     1613
     1614
     1615
     1616    def test_bedslope_problem_second_order_more_steps(self):
     1617
     1618        from mesh_factory import rectangular
     1619        from shallow_water import Domain, Reflective_boundary, Constant_height
     1620        from Numeric import array
     1621
     1622        #Create basic mesh
     1623        points, vertices, boundary = rectangular(6, 6)
     1624
     1625        #Create shallow water domain
     1626        domain = Domain(points, vertices, boundary)
     1627        domain.smooth = False
     1628        domain.default_order=2
     1629
     1630        #Bed-slope and friction at vertices (and interpolated elsewhere)
     1631        def x_slope(x, y):
     1632            return -x/3
     1633
     1634        domain.set_quantity('elevation', x_slope)
     1635
     1636        # Boundary conditions
     1637        Br = Reflective_boundary(domain)
     1638        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
     1639
     1640        #Initial condition
     1641        domain.set_quantity('level', Constant_height(x_slope, 0.05))
     1642        domain.check_integrity()
     1643
     1644        assert allclose(domain.quantities['level'].centroid_values,
     1645                        [0.01296296, 0.03148148, 0.01296296,
     1646                        0.03148148, 0.01296296, 0.03148148,
     1647                        0.01296296, 0.03148148, 0.01296296,
     1648                        0.03148148, 0.01296296, 0.03148148,
     1649                        -0.04259259, -0.02407407, -0.04259259,
     1650                        -0.02407407, -0.04259259, -0.02407407,
     1651                        -0.04259259, -0.02407407, -0.04259259,
     1652                        -0.02407407, -0.04259259, -0.02407407,
     1653                        -0.09814815, -0.07962963, -0.09814815,
     1654                        -0.07962963, -0.09814815, -0.07962963,
     1655                        -0.09814815, -0.07962963, -0.09814815,
     1656                        -0.07962963, -0.09814815, -0.07962963,
     1657                        -0.1537037 , -0.13518519, -0.1537037,
     1658                        -0.13518519, -0.1537037, -0.13518519,
     1659                        -0.1537037 , -0.13518519, -0.1537037,
     1660                        -0.13518519, -0.1537037, -0.13518519,
     1661                        -0.20925926, -0.19074074, -0.20925926,
     1662                        -0.19074074, -0.20925926, -0.19074074,
     1663                        -0.20925926, -0.19074074, -0.20925926,
     1664                        -0.19074074, -0.20925926, -0.19074074,
     1665                        -0.26481481, -0.2462963, -0.26481481,
     1666                        -0.2462963, -0.26481481, -0.2462963,
     1667                        -0.26481481, -0.2462963, -0.26481481,
     1668                        -0.2462963, -0.26481481, -0.2462963])
     1669
     1670
     1671        #print domain.quantities['level'].extrapolate_second_order()
     1672        #domain.distribute_to_vertices_and_edges()
     1673        #print domain.quantities['level'].vertex_values[:,0]       
     1674       
     1675        #Evolution
     1676        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
     1677            pass
     1678
     1679       
     1680        assert allclose(domain.quantities['level'].centroid_values,
     1681     [-0.02907028, -0.01475478, -0.02973417, -0.01447186, -0.02932665, -0.01428285,
     1682      -0.02901975, -0.0141361,  -0.02898816, -0.01418135, -0.02961409, -0.01403487,
     1683      -0.07597998, -0.06252591, -0.07664854, -0.06312532, -0.07638287, -0.06265139,
     1684      -0.07571145, -0.06235231, -0.0756817,  -0.06245309, -0.07652292, -0.06289946,
     1685      -0.12367464, -0.11088981, -0.12237277, -0.11115338, -0.1218934,  -0.1107174,
     1686      -0.12081485, -0.11000491, -0.12038451, -0.11010335, -0.12102113, -0.11012105,
     1687      -0.16909116, -0.15831543, -0.16730214, -0.15786249, -0.1665493,  -0.15697919,
     1688      -0.16496618, -0.15559852, -0.16338679, -0.15509088, -0.16364092, -0.15424423,
     1689      -0.18771107, -0.19903904, -0.18903759, -0.19858437, -0.18701552, -0.19697797,
     1690      -0.1833593,  -0.19505871, -0.1818806,  -0.19418042, -0.18586159, -0.19576946,
     1691      -0.13986873, -0.14170053, -0.14132188, -0.14560674, -0.14095617, -0.14373292,
     1692      -0.13785933, -0.14033364, -0.13592955, -0.13936356, -0.13596008, -0.14216296])
     1693
     1694        assert allclose(domain.quantities['xmomentum'].centroid_values,
     1695     [ 0.00831121,  0.00317948,  0.00731797,  0.00334939,  0.00764717,  0.00348053,
     1696       0.00788729,  0.00356522,  0.00780649,  0.00341919,  0.00693525,  0.00310375,
     1697       0.02166196,  0.01421475,  0.02017737,  0.01316839,  0.02037015,  0.01368659,
     1698       0.02106,     0.01399161,  0.02074514,  0.01354935,  0.01887407,  0.0123113,
     1699       0.03775083,  0.02855197,  0.03689337,  0.02759782,  0.03732848,  0.02812072,
     1700       0.03872545,  0.02913348,  0.03880939,  0.02803804,  0.03546499,  0.0260039,
     1701       0.0632131,   0.04730634,  0.0576324,   0.04592336,  0.05790921,  0.04690514,
     1702       0.05986467,  0.04871165,  0.06170068,  0.04811572,  0.05657041,  0.04416292,
     1703       0.08489642,  0.07188097,  0.07835261,  0.06843406,  0.07986412,  0.0698247,
     1704       0.08201071,  0.07216756,  0.08378418,  0.07273624,  0.080399,    0.06645841,
     1705       0.01631548,  0.04691608,  0.0206632,   0.044409,    0.02115518,  0.04560305,
     1706       0.02160608,  0.04663725,  0.02174734,  0.04795559,  0.02281427,  0.05667111])
     1707
     1708       
     1709        assert allclose(domain.quantities['ymomentum'].centroid_values,
     1710     [ 1.45876601e-004, -3.24627393e-004, -1.57572719e-004, -2.92790187e-004,
     1711      -9.90988382e-005, -3.06677335e-004, -1.62493106e-004, -3.71310004e-004,
     1712      -1.99445058e-004, -3.28493467e-004,  6.68217349e-005, -8.42042805e-006,
     1713       5.05093371e-004, -1.42842214e-004, -6.81454718e-005, -5.02084057e-004,
     1714      -8.50583861e-005, -4.65443981e-004, -1.96406564e-004, -5.88889562e-004,
     1715      -2.70160173e-004, -5.35485454e-004,  2.60780997e-004,  3.12145471e-005,
     1716       5.16189608e-004,  1.07069062e-004,  9.29989252e-005, -3.71211119e-004,
     1717       1.16350246e-004, -3.82407830e-004, -1.62077969e-004, -6.30906636e-004,
     1718      -4.74025708e-004, -6.94463009e-004,  6.15092843e-005,  2.22106820e-004,
     1719      -6.29589294e-004,  2.43611937e-004, -5.88125094e-004, -6.94293192e-005,
     1720      -4.17914641e-004,  6.64609019e-005, -7.68334577e-004, -3.40232101e-004,
     1721      -1.67424308e-003, -7.39485066e-004, -1.59966988e-003,  5.68262838e-005,
     1722      -1.48470633e-003, -1.84554882e-003, -2.27200099e-003, -1.67506848e-003,
     1723      -1.95610258e-003, -1.47638801e-003, -1.73779477e-003, -1.85498791e-003,
     1724      -2.01357843e-003, -2.17675471e-003, -1.65783870e-003, -1.15818681e-003,
     1725      -1.18663036e-003, -2.94229849e-003, -3.59309018e-003, -5.13496584e-003,
     1726      -6.17359400e-003, -5.98761937e-003, -6.00540116e-003, -5.01121966e-003,
     1727      -4.50964850e-003, -3.06319963e-003,  6.08950810e-004, -4.79537921e-004])
     1728       
     1729
     1730       
     1731
     1732    def test_temp_play(self):
     1733
     1734        from mesh_factory import rectangular
     1735        from shallow_water import Domain, Reflective_boundary, Constant_height
     1736        from Numeric import array
     1737
     1738        #Create basic mesh
     1739        points, vertices, boundary = rectangular(5, 5)
     1740
     1741        #Create shallow water domain
     1742        domain = Domain(points, vertices, boundary)
     1743        domain.smooth = False
     1744        domain.default_order=2
     1745
     1746        #Bed-slope and friction at vertices (and interpolated elsewhere)
     1747        def x_slope(x, y):
     1748            return -x/3
     1749
     1750        domain.set_quantity('elevation', x_slope)
     1751
     1752        # Boundary conditions
     1753        Br = Reflective_boundary(domain)
     1754        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
     1755
     1756        #Initial condition
     1757        domain.set_quantity('level', Constant_height(x_slope, 0.05))
     1758        domain.check_integrity()
     1759
     1760        #Evolution
     1761        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1):
     1762            pass
     1763
     1764        assert allclose(domain.quantities['level'].centroid_values[:4],
     1765                        [0.00206836, 0.01296714, 0.00363415, 0.01438924])
     1766        assert allclose(domain.quantities['xmomentum'].centroid_values[:4],                     
     1767                        [0.01360154, 0.00671133, 0.01264578, 0.00648503])
     1768        assert allclose(domain.quantities['ymomentum'].centroid_values[:4],
     1769                        [-1.19201077e-003, -7.23647546e-004,
     1770                        -6.39083123e-005, 6.29815168e-005])
     1771       
     1772       
     1773       
    13911774       
    13921775   
Note: See TracChangeset for help on using the changeset viewer.