Changeset 7731
- Timestamp:
- May 18, 2010, 2:54:05 PM (15 years ago)
- Location:
- anuga_core/source/anuga/shallow_water
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/shallow_water/__init__.py
r7519 r7731 7 7 8 8 # Make selected classes available directly 9 from shallow_water_domain import Domain,\ 10 Transmissive_boundary, Reflective_boundary,\ 11 Dirichlet_boundary, Time_boundary, File_boundary,\ 9 from boundaries import Reflective_boundary,\ 12 10 Transmissive_momentum_set_stage_boundary,\ 13 11 Dirichlet_discharge_boundary,\ 14 12 Field_boundary,\ 15 13 Transmissive_stage_zero_momentum_boundary,\ 16 Transmissive_n_momentum_zero_t_momentum_set_stage_boundary,\ 17 AWI_boundary 14 Transmissive_n_momentum_zero_t_momentum_set_stage_boundary 15 16 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\ 17 import Transmissive_boundary, Dirichlet_boundary, \ 18 Time_boundary, File_boundary, AWI_boundary 19 20 from shallow_water_domain import Domain 18 21 19 22 #from shallow_water_balanced_domain import Swb_domain 20 23 21 24 22 # FIXME (Ole): Deprecate23 from shallow_water_domain import Transmissive_Momentum_Set_Stage_boundary24 from shallow_water_domain import Dirichlet_Discharge_boundary25 25 -
anuga_core/source/anuga/shallow_water/shallow_water_domain.py
r7711 r7731 1213 1213 1214 1214 1215 ################################################################################1216 # Boundary conditions - specific to the shallow water wave equation1217 ################################################################################1218 1219 ##1220 # @brief Class for a reflective boundary.1221 # @note Inherits from Boundary.1222 class Reflective_boundary(Boundary):1223 """Reflective boundary returns same conserved quantities as1224 those present in its neighbour volume but reflected.1225 1226 This class is specific to the shallow water equation as it1227 works with the momentum quantities assumed to be the second1228 and third conserved quantities.1229 """1230 1231 ##1232 # @brief Instantiate a Reflective_boundary.1233 # @param domain1234 def __init__(self, domain=None):1235 Boundary.__init__(self)1236 1237 if domain is None:1238 msg = 'Domain must be specified for reflective boundary'1239 raise Exception, msg1240 1241 # Handy shorthands1242 self.stage = domain.quantities['stage'].edge_values1243 self.xmom = domain.quantities['xmomentum'].edge_values1244 self.ymom = domain.quantities['ymomentum'].edge_values1245 self.normals = domain.normals1246 1247 self.conserved_quantities = num.zeros(3, num.float)1248 1249 ##1250 # @brief Return a representation of this instance.1251 def __repr__(self):1252 return 'Reflective_boundary'1253 1254 ##1255 # @brief Calculate reflections (reverse outward momentum).1256 # @param vol_id1257 # @param edge_id1258 def evaluate(self, vol_id, edge_id):1259 """Reflective boundaries reverses the outward momentum1260 of the volume they serve.1261 """1262 1263 q = self.conserved_quantities1264 q[0] = self.stage[vol_id, edge_id]1265 q[1] = self.xmom[vol_id, edge_id]1266 q[2] = self.ymom[vol_id, edge_id]1267 1268 normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]1269 1270 r = rotate(q, normal, direction = 1)1271 r[1] = -r[1]1272 q = rotate(r, normal, direction = -1)1273 1274 return q1275 1276 1277 ##1278 # @brief Class for a transmissive boundary.1279 # @note Inherits from Boundary.1280 class Transmissive_momentum_set_stage_boundary(Boundary):1281 """Returns same momentum conserved quantities as1282 those present in its neighbour volume.1283 Sets stage by specifying a function f of time which may either be a1284 vector function or a scalar function1285 1286 Example:1287 1288 def waveform(t):1289 return sea_level + normalized_amplitude/cosh(t-25)**21290 1291 Bts = Transmissive_momentum_set_stage_boundary(domain, waveform)1292 1293 Underlying domain must be specified when boundary is instantiated1294 """1295 1296 ##1297 # @brief Instantiate a Transmissive_momentum_set_stage_boundary.1298 # @param domain1299 # @param function1300 def __init__(self, domain=None, function=None):1301 Boundary.__init__(self)1302 1303 if domain is None:1304 msg = 'Domain must be specified for this type boundary'1305 raise Exception, msg1306 1307 if function is None:1308 msg = 'Function must be specified for this type boundary'1309 raise Exception, msg1310 1311 self.domain = domain1312 self.function = function1313 1314 ##1315 # @brief Return a representation of this instance.1316 def __repr__(self):1317 return 'Transmissive_momentum_set_stage_boundary(%s)' %self.domain1318 1319 ##1320 # @brief Calculate transmissive results.1321 # @param vol_id1322 # @param edge_id1323 def evaluate(self, vol_id, edge_id):1324 """Transmissive momentum set stage boundaries return the edge momentum1325 values of the volume they serve.1326 """1327 1328 q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)1329 1330 t = self.domain.get_time()1331 1332 if hasattr(self.function, 'time'):1333 # Roll boundary over if time exceeds1334 while t > self.function.time[-1]:1335 msg = 'WARNING: domain time %.2f has exceeded' % t1336 msg += 'time provided in '1337 msg += 'transmissive_momentum_set_stage_boundary object.\n'1338 msg += 'I will continue, reusing the object from t==0'1339 log.critical(msg)1340 t -= self.function.time[-1]1341 1342 value = self.function(t)1343 try:1344 x = float(value)1345 except:1346 x = float(value[0])1347 1348 q[0] = x1349 1350 return q1351 1352 # FIXME: Consider this (taken from File_boundary) to allow1353 # spatial variation1354 # if vol_id is not None and edge_id is not None:1355 # i = self.boundary_indices[ vol_id, edge_id ]1356 # return self.F(t, point_id = i)1357 # else:1358 # return self.F(t)1359 1360 1361 ##1362 # @brief Deprecated boundary class.1363 class Transmissive_Momentum_Set_Stage_boundary(Transmissive_momentum_set_stage_boundary):1364 pass1365 1366 1367 ##1368 # @brief Class for a transmissive boundary.1369 # @note Inherits from Boundary.1370 class Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(Boundary):1371 """Returns the same normal momentum as that1372 present in neighbour volume edge. Zero out the tangential momentum.1373 Sets stage by specifying a function f of time which may either be a1374 vector function or a scalar function1375 1376 Example:1377 1378 def waveform(t):1379 return sea_level + normalized_amplitude/cosh(t-25)**21380 1381 Bts = Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(domain, waveform)1382 1383 Underlying domain must be specified when boundary is instantiated1384 """1385 1386 ##1387 # @brief Instantiate a Transmissive_n_momentum_zero_t_momentum_set_stage_boundary.1388 # @param domain1389 # @param function1390 def __init__(self, domain=None, function=None):1391 Boundary.__init__(self)1392 1393 if domain is None:1394 msg = 'Domain must be specified for this type boundary'1395 raise Exception, msg1396 1397 if function is None:1398 msg = 'Function must be specified for this type boundary'1399 raise Exception, msg1400 1401 self.domain = domain1402 self.function = function1403 1404 ##1405 # @brief Return a representation of this instance.1406 def __repr__(self):1407 return 'Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(%s)' %self.domain1408 1409 ##1410 # @brief Calculate transmissive results.1411 # @param vol_id1412 # @param edge_id1413 def evaluate(self, vol_id, edge_id):1414 """Transmissive_n_momentum_zero_t_momentum_set_stage_boundary1415 return the edge momentum values of the volume they serve.1416 """1417 1418 q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)1419 1420 normal = self.domain.get_normal(vol_id, edge_id)1421 1422 1423 t = self.domain.get_time()1424 1425 if hasattr(self.function, 'time'):1426 # Roll boundary over if time exceeds1427 while t > self.function.time[-1]:1428 msg = 'WARNING: domain time %.2f has exceeded' % t1429 msg += 'time provided in '1430 msg += 'transmissive_momentum_set_stage_boundary object.\n'1431 msg += 'I will continue, reusing the object from t==0'1432 log.critical(msg)1433 t -= self.function.time[-1]1434 1435 value = self.function(t)1436 try:1437 x = float(value)1438 except:1439 x = float(value[0])1440 1441 ## import math1442 ## if vol_id == 9433:1443 ## print 'vol_id = ',vol_id, ' edge_id = ',edge_id, 'q = ', q, ' x = ',x1444 ## print 'normal = ', normal1445 ## print 'n . p = ', (normal[0]*q[1] + normal[1]*q[2])1446 ## print 't . p = ', (normal[1]*q[1] - normal[0]*q[2])1447 1448 1449 q[0] = x1450 ndotq = (normal[0]*q[1] + normal[1]*q[2])1451 q[1] = normal[0]*ndotq1452 q[2] = normal[1]*ndotq1453 1454 1455 return q1456 1457 ##1458 # @brief A transmissive boundary, momentum set to zero.1459 # @note Inherits from Bouondary.1460 class Transmissive_stage_zero_momentum_boundary(Boundary):1461 """Return same stage as those present in its neighbour volume.1462 Set momentum to zero.1463 1464 Underlying domain must be specified when boundary is instantiated1465 """1466 1467 ##1468 # @brief Instantiate a Transmissive (zero momentum) boundary.1469 # @param domain1470 def __init__(self, domain=None):1471 Boundary.__init__(self)1472 1473 if domain is None:1474 msg = ('Domain must be specified for '1475 'Transmissive_stage_zero_momentum boundary')1476 raise Exception, msg1477 1478 self.domain = domain1479 1480 ##1481 # @brief Return a representation of this instance.1482 def __repr__(self):1483 return 'Transmissive_stage_zero_momentum_boundary(%s)' % self.domain1484 1485 ##1486 # @brief Calculate transmissive (zero momentum) results.1487 # @param vol_id1488 # @param edge_id1489 def evaluate(self, vol_id, edge_id):1490 """Transmissive boundaries return the edge values1491 of the volume they serve.1492 """1493 1494 q = self.domain.get_conserved_quantities(vol_id, edge=edge_id)1495 1496 q[1] = q[2] = 0.01497 return q1498 1499 1500 ##1501 # @brief Class for a Dirichlet discharge boundary.1502 # @note Inherits from Boundary.1503 class Dirichlet_discharge_boundary(Boundary):1504 """1505 Sets stage (stage0)1506 Sets momentum (wh0) in the inward normal direction.1507 1508 Underlying domain must be specified when boundary is instantiated1509 """1510 1511 ##1512 # @brief Instantiate a Dirichlet discharge boundary.1513 # @param domain1514 # @param stage01515 # @param wh01516 def __init__(self, domain=None, stage0=None, wh0=None):1517 Boundary.__init__(self)1518 1519 if domain is None:1520 msg = 'Domain must be specified for this type of boundary'1521 raise Exception, msg1522 1523 if stage0 is None:1524 raise Exception, 'Stage must be specified for this type of boundary'1525 1526 if wh0 is None:1527 wh0 = 0.01528 1529 self.domain = domain1530 self.stage0 = stage01531 self.wh0 = wh01532 1533 ##1534 # @brief Return a representation of this instance.1535 def __repr__(self):1536 return 'Dirichlet_Discharge_boundary(%s)' % self.domain1537 1538 ##1539 # @brief Calculate Dirichlet discharge boundary results.1540 # @param vol_id1541 # @param edge_id1542 def evaluate(self, vol_id, edge_id):1543 """Set discharge in the (inward) normal direction"""1544 1545 normal = self.domain.get_normal(vol_id,edge_id)1546 q = [self.stage0, -self.wh0*normal[0], -self.wh0*normal[1]]1547 return q1548 1549 # FIXME: Consider this (taken from File_boundary) to allow1550 # spatial variation1551 # if vol_id is not None and edge_id is not None:1552 # i = self.boundary_indices[ vol_id, edge_id ]1553 # return self.F(t, point_id = i)1554 # else:1555 # return self.F(t)1556 1557 1558 # Backward compatibility1559 # FIXME(Ole): Deprecate1560 ##1561 # @brief Deprecated1562 class Dirichlet_Discharge_boundary(Dirichlet_discharge_boundary):1563 pass1564 1565 1566 class Inflow_boundary(Boundary):1567 """Apply given flow in m^3/s to boundary segment.1568 Depth and momentum is derived using Manning's formula.1569 1570 Underlying domain must be specified when boundary is instantiated1571 """1572 1573 # FIXME (Ole): This is work in progress and definitely not finished.1574 # The associated test has been disabled1575 1576 def __init__(self, domain=None, rate=0.0):1577 Boundary.__init__(self)1578 1579 if domain is None:1580 msg = 'Domain must be specified for '1581 msg += 'Inflow boundary'1582 raise Exception, msg1583 1584 self.domain = domain1585 1586 # FIXME(Ole): Allow rate to be time dependent as well1587 self.rate = rate1588 self.tag = None # Placeholder for tag associated with this object.1589 1590 def __repr__(self):1591 return 'Inflow_boundary(%s)' %self.domain1592 1593 def evaluate(self, vol_id, edge_id):1594 """Apply inflow rate at each edge of this boundary1595 """1596 1597 # First find all segments having the same tag is vol_id, edge_id1598 # This will be done the first time evaluate is called.1599 if self.tag is None:1600 boundary = self.domain.boundary1601 self.tag = boundary[(vol_id, edge_id)]1602 1603 # Find total length of boundary with this tag1604 length = 0.01605 for v_id, e_id in boundary:1606 if self.tag == boundary[(v_id, e_id)]:1607 length += self.domain.mesh.get_edgelength(v_id, e_id)1608 1609 self.length = length1610 self.average_momentum = self.rate/length1611 1612 1613 # Average momentum has now been established across this boundary1614 # Compute momentum in the inward normal direction1615 1616 inward_normal = -self.domain.mesh.get_normal(vol_id, edge_id)1617 xmomentum, ymomentum = self.average_momentum * inward_normal1618 1619 # Compute depth based on Manning's formula v = 1/n h^{2/3} sqrt(S)1620 # Where v is velocity, n is manning's coefficient, h is depth and S is the slope into the domain.1621 # Let mu be the momentum (vh), then this equation becomes: mu = 1/n h^{5/3} sqrt(S)1622 # from which we can isolate depth to get1623 # h = (mu n/sqrt(S) )^{3/5}1624 1625 slope = 0 # get gradient for this triangle dot normal1626 1627 # get manning coef from this triangle1628 friction = self.domain.get_quantity('friction').get_values(location='edges',1629 indices=[vol_id])[0]1630 mannings_n = friction[edge_id]1631 1632 if slope > epsilon and mannings_n > epsilon:1633 depth = pow(self.average_momentum * mannings_n/math.sqrt(slope), 3.0/5)1634 else:1635 depth = 1.01636 1637 # Elevation on this edge1638 1639 z = self.domain.get_quantity('elevation').get_values(location='edges',1640 indices=[vol_id])[0]1641 elevation = z[edge_id]1642 1643 # Assign conserved quantities and return1644 q = num.array([elevation + depth, xmomentum, ymomentum], num.Float)1645 return q1646 1647 1648 1649 1650 1651 1652 class Field_boundary(Boundary):1653 """Set boundary from given field represented in an sww file containing1654 values for stage, xmomentum and ymomentum.1655 1656 Optionally, the user can specify mean_stage to offset the stage provided1657 in the sww file.1658 1659 This function is a thin wrapper around the generic File_boundary. The1660 difference between the file_boundary and field_boundary is only that the1661 field_boundary will allow you to change the level of the stage height when1662 you read in the boundary condition. This is very useful when running1663 different tide heights in the same area as you need only to convert one1664 boundary condition to a SWW file, ideally for tide height of 0 m1665 (saving disk space). Then you can use field_boundary to read this SWW file1666 and change the stage height (tide) on the fly depending on the scenario.1667 """1668 1669 ##1670 # @brief Construct an instance of a 'field' boundary.1671 # @param filename Name of SWW file containing stage, x and ymomentum.1672 # @param domain Shallow water domain for which the boundary applies.1673 # @param mean_stage Mean water level added to stage derived from SWW file.1674 # @param time_thinning Time step thinning factor.1675 # @param time_limit1676 # @param boundary_polygon1677 # @param default_boundary None or an instance of Boundary.1678 # @param use_cache True if caching is to be used.1679 # @param verbose True if this method is to be verbose.1680 def __init__(self,1681 filename,1682 domain,1683 mean_stage=0.0,1684 time_thinning=1,1685 time_limit=None,1686 boundary_polygon=None,1687 default_boundary=None,1688 use_cache=False,1689 verbose=False):1690 """Constructor1691 1692 filename: Name of sww file1693 domain: pointer to shallow water domain for which the boundary applies1694 mean_stage: The mean water level which will be added to stage derived1695 from the boundary condition1696 time_thinning: Will set how many time steps from the sww file read in1697 will be interpolated to the boundary. For example if1698 the sww file has 1 second time steps and is 24 hours1699 in length it has 86400 time steps. If you set1700 time_thinning to 1 it will read all these steps.1701 If you set it to 100 it will read every 100th step eg1702 only 864 step. This parameter is very useful to increase1703 the speed of a model run that you are setting up1704 and testing.1705 1706 default_boundary: Must be either None or an instance of a1707 class descending from class Boundary.1708 This will be used in case model time exceeds1709 that available in the underlying data.1710 1711 Note that mean_stage will also be added to this.1712 1713 use_cache:1714 verbose:1715 """1716 1717 # Create generic file_boundary object1718 self.file_boundary = File_boundary(filename,1719 domain,1720 time_thinning=time_thinning,1721 time_limit=time_limit,1722 boundary_polygon=boundary_polygon,1723 default_boundary=default_boundary,1724 use_cache=use_cache,1725 verbose=verbose)1726 1727 # Record information from File_boundary1728 self.F = self.file_boundary.F1729 self.domain = self.file_boundary.domain1730 1731 # Record mean stage1732 self.mean_stage = mean_stage1733 1734 ##1735 # @note Generate a string representation of this instance.1736 def __repr__(self):1737 return 'Field boundary'1738 1739 ##1740 # @brief Calculate 'field' boundary results.1741 # @param vol_id1742 # @param edge_id1743 def evaluate(self, vol_id=None, edge_id=None):1744 """Return linearly interpolated values based on domain.time1745 1746 vol_id and edge_id are ignored1747 """1748 1749 # Evaluate file boundary1750 q = self.file_boundary.evaluate(vol_id, edge_id)1751 1752 # Adjust stage1753 for j, name in enumerate(self.domain.conserved_quantities):1754 if name == 'stage':1755 q[j] += self.mean_stage1756 return q1757 1758 1215 1759 1216 ################################################################################ … … 2826 2283 if compile.can_use_C_extension('shallow_water_ext.c'): 2827 2284 # Underlying C implementations can be accessed 2828 from shallow_water_ext import rotate,assign_windfield_values2285 from shallow_water_ext import assign_windfield_values 2829 2286 else: 2830 2287 msg = 'C implementations could not be accessed by %s.\n ' % __file__ -
anuga_core/source/anuga/shallow_water/test_read_sww.py
r7562 r7731 37 37 from anuga.interface import rectangular_cross 38 38 from anuga.interface import Domain 39 from anuga.interface import Reflective_boundary 40 from anuga.interface import Dirichlet_boundary 39 from boundaries import Reflective_boundary 40 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\ 41 import Dirichlet_boundary 41 42 from anuga.interface import Time_boundary 42 43 -
anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py
r7718 r7731 8 8 from anuga.config import g, epsilon 9 9 from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a 10 from anuga.utilities.numerical_tools import mean 10 from anuga.utilities.numerical_tools import mean, ensure_numeric 11 11 from anuga.geometry.polygon import is_inside_polygon 12 12 from anuga.coordinate_transforms.geo_reference import Geo_reference 13 from anuga.abstract_2d_finite_volumes.quantity import Quantity14 13 from anuga.geospatial_data.geospatial_data import Geospatial_data 15 14 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross 15 from anuga.abstract_2d_finite_volumes.quantity import Quantity 16 16 17 17 from anuga.utilities.system_tools import get_pathname_from_package 18 from shallow_water_domain import * 18 19 from anuga.shallow_water import Domain 20 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions \ 21 import Dirichlet_boundary 22 from anuga.shallow_water.shallow_water_domain import Rainfall, Wind_stress 23 from anuga.shallow_water.shallow_water_domain import Inflow, Cross_section 24 from anuga.shallow_water.data_manager import get_flow_through_cross_section 25 26 # boundary functions 27 from anuga.shallow_water.boundaries import Reflective_boundary, \ 28 Field_boundary, Transmissive_momentum_set_stage_boundary, \ 29 Transmissive_stage_zero_momentum_boundary 30 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\ 31 import Transmissive_boundary, Dirichlet_boundary, \ 32 Time_boundary, File_boundary, AWI_boundary 19 33 20 34 import numpy as num … … 22 36 # Get gateway to C implementation of flux function for direct testing 23 37 from shallow_water_ext import flux_function_central as flux_function 38 from shallow_water_ext import rotate 24 39 25 40 … … 6601 6616 import rectangular_cross 6602 6617 from anuga.shallow_water import Domain 6603 from anuga.shallow_water.shallow_water_domain import Reflective_boundary6604 from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary6605 6618 from anuga.shallow_water.shallow_water_domain import Inflow 6606 6619 from anuga.shallow_water.data_manager \ … … 6694 6707 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross 6695 6708 from anuga.shallow_water import Domain 6696 from anuga.shallow_water.shallow_water_domain import Reflective_boundary6697 from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary6698 6709 from anuga.shallow_water.shallow_water_domain import Inflow 6699 6710 from anuga.shallow_water.data_manager import get_flow_through_cross_section … … 6790 6801 verbose = False 6791 6802 6792 6793 #---------------------------------------------------------------------6794 # Import necessary modules6795 #---------------------------------------------------------------------6796 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross6797 from anuga.shallow_water import Domain6798 from anuga.shallow_water.shallow_water_domain import Reflective_boundary6799 from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary6800 from anuga.shallow_water.shallow_water_domain import Rainfall6801 from anuga.shallow_water.data_manager import get_flow_through_cross_section6802 6803 6803 6804 #---------------------------------------------------------------------- … … 7125 7126 verbose = False 7126 7127 7127 7128 #----------------------------------------------------------------------7129 # Import necessary modules7130 #----------------------------------------------------------------------7131 7132 from anuga.abstract_2d_finite_volumes.mesh_factory \7133 import rectangular_cross7134 from anuga.shallow_water import Domain7135 from anuga.shallow_water.shallow_water_domain import Reflective_boundary7136 from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary7137 from anuga.shallow_water.shallow_water_domain import Inflow7138 from anuga.shallow_water.data_manager \7139 import get_flow_through_cross_section7140 from anuga.abstract_2d_finite_volumes.util \7141 import sww2csv_gauges, csv2timeseries_graphs7142 7143 7128 #---------------------------------------------------------------------- 7144 7129 # Setup computational domain
Note: See TracChangeset
for help on using the changeset viewer.