Changeset 229
- Timestamp:
- Aug 26, 2004, 10:14:19 PM (21 years ago)
- Location:
- inundation/ga/storm_surge/pyvolution
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/ga/storm_surge/pyvolution/domain.py
r218 r229 290 290 #Update boundary values 291 291 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] 292 302 293 303 #Compute all fluxes and timestep suitable for all volumes 294 304 self.compute_fluxes() 305 ##print 306 ##for name in self.conserved_quantities: 307 ## Q = self.quantities[name] 308 ## print 'EU:', Q.explicit_update[:4] 295 309 296 310 #Update timestep to fit yieldstep and finaltime … … 300 314 self.update_conserved_quantities() 301 315 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 302 326 #Update vertex and edge values 303 327 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 304 336 305 337 #Update time … … 310 342 self.number_of_first_order_steps += 1 311 343 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 312 351 #Yield results 313 352 if finaltime is not None and abs(self.time - finaltime) < epsilon: -
inundation/ga/storm_surge/pyvolution/quantity.py
r221 r229 78 78 v2 = self.vertex_values[i, 2] 79 79 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' 82 83 self.edge_values[i, 0] = 0.5*(v1 + v2) 83 84 self.edge_values[i, 1] = 0.5*(v0 + v2) -
inundation/ga/storm_surge/pyvolution/shallow_water.py
r221 r229 273 273 def distribute_to_vertices_and_edges(domain): 274 274 275 # print 'Calling distrib within SW'275 ##print 'Calling distrib within SW' 276 276 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) 277 280 protect_against_infinitesimal_heights_centroid(domain) 278 281 extrapolate_first_order(domain) 279 282 elif domain.order == 2: 280 protect_against_negative_heights_centroid(domain)283 #protect_against_negative_heights_centroid(domain) 281 284 extrapolate_second_order(domain) 282 285 else: … … 323 326 """ 324 327 325 #FIXME: USed only in second order328 326 329 #Water levels at centroids 327 330 wc = domain.quantities['level'].centroid_values … … 428 431 #Update water levels 429 432 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 434 439 435 440 def extrapolate_second_order(domain): … … 462 467 Q = domain.quantities[name] 463 468 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] 464 479 Q.limit() 465 480 466 481 482 #print 'y1', Q.vertex_values[1,:] #OK 483 467 484 #Water levels at centroids 468 485 wc = domain.quantities['level'].centroid_values … … 507 524 #first order scheme and alpha==1 means using the limited 508 525 #second order scheme for the stage w 509 #If hmin < 0 then alpha=0 rever rting to first order.526 #If hmin < 0 then alpha=0 reverting to first order. 510 527 511 528 if z_range > 0.0: … … 513 530 else: 514 531 alpha = 1.0 515 532 533 ##if k==1: print 'alpha', alpha #OK 534 516 535 #Update water levels 517 536 # (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 538 560 #Finally, protect against infinitesimal heights and high speeds 539 561 #Water depths at vertices -
inundation/ga/storm_surge/pyvolution/test_shallow_water.py
r221 r229 859 859 assert allclose(L[1], [4.1, 16.1, 20.1]) 860 860 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 862 912 863 913 def test_second_order_flat_bed_onestep(self): … … 1297 1347 -0.14183568, -0.13623283, -0.14529803]) 1298 1348 1299 def test_bedslope_problem_second_order (self):1349 def test_bedslope_problem_second_order_one_step(self): 1300 1350 1301 1351 from mesh_factory import rectangular … … 1357 1407 1358 1408 #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): 1360 1410 #domain.write_time() 1361 1411 pass 1362 1412 1363 1413 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 1367 1501 assert allclose(domain.min_timestep, 0.0376895634803) 1368 1502 assert allclose(domain.max_timestep, 0.0415635655309) 1369 1503 1370 1504 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 1391 1774 1392 1775
Note: See TracChangeset
for help on using the changeset viewer.