Changeset 7176 for branches/numpy/anuga/shallow_water
- Timestamp:
- Jun 10, 2009, 5:28:35 PM (16 years ago)
- Location:
- branches/numpy/anuga/shallow_water
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/numpy/anuga/shallow_water/data_manager.py
r7035 r7176 566 566 # Define a zero vector of same size and type as A 567 567 # for use with momenta 568 null = num.zeros(num.size(A), A.dtype.char) #??#568 null = num.zeros(num.size(A), A.dtype.char) 569 569 570 570 # Get xmomentum where depth exceeds minimum_storable_height … … 1433 1433 for name in infile.variables: 1434 1434 var = infile.variables[name] 1435 outfile.createVariable(name, var.dtype.char, var.dimensions) #??#1435 outfile.createVariable(name, var.dtype.char, var.dimensions) 1436 1436 1437 1437 # Copy the static variables … … 2193 2193 swwfile = basename_in + '.sww' 2194 2194 demfile = basename_out + '.' + format 2195 # Note the use of a .ers extension is optional (write_ermapper_grid will2196 # deal with either option2197 2195 2198 2196 # Read sww file … … 2697 2695 # Interpolate using quantity values 2698 2696 if verbose: print 'Interpolating' 2699 interpolated_values = interp.interpolate(q, data_points).flatten() #????#2697 interpolated_values = interp.interpolate(q, data_points).flatten() 2700 2698 2701 2699 if verbose: -
branches/numpy/anuga/shallow_water/shallow_water_domain.py
r7035 r7176 1462 1462 1463 1463 1464 ## 1465 # @brief A class for a 'field' boundary. 1466 # @note Inherits from Boundary. 1464 class Inflow_boundary(Boundary): 1465 """Apply given flow in m^3/s to boundary segment. 1466 Depth and momentum is derived using Manning's formula. 1467 1468 Underlying domain must be specified when boundary is instantiated 1469 """ 1470 1471 # FIXME (Ole): This is work in progress and definitely not finished. 1472 # The associated test has been disabled 1473 1474 def __init__(self, domain=None, rate=0.0): 1475 Boundary.__init__(self) 1476 1477 if domain is None: 1478 msg = 'Domain must be specified for ' 1479 msg += 'Inflow boundary' 1480 raise Exception, msg 1481 1482 self.domain = domain 1483 1484 # FIXME(Ole): Allow rate to be time dependent as well 1485 self.rate = rate 1486 self.tag = None # Placeholder for tag associated with this object. 1487 1488 def __repr__(self): 1489 return 'Inflow_boundary(%s)' %self.domain 1490 1491 def evaluate(self, vol_id, edge_id): 1492 """Apply inflow rate at each edge of this boundary 1493 """ 1494 1495 # First find all segments having the same tag is vol_id, edge_id 1496 # This will be done the first time evaluate is called. 1497 if self.tag is None: 1498 boundary = self.domain.boundary 1499 self.tag = boundary[(vol_id, edge_id)] 1500 1501 # Find total length of boundary with this tag 1502 length = 0.0 1503 for v_id, e_id in boundary: 1504 if self.tag == boundary[(v_id, e_id)]: 1505 length += self.domain.mesh.get_edgelength(v_id, e_id) 1506 1507 self.length = length 1508 self.average_momentum = self.rate/length 1509 1510 1511 # Average momentum has now been established across this boundary 1512 # Compute momentum in the inward normal direction 1513 1514 inward_normal = -self.domain.mesh.get_normal(vol_id, edge_id) 1515 xmomentum, ymomentum = self.average_momentum * inward_normal 1516 1517 # Compute depth based on Manning's formula v = 1/n h^{2/3} sqrt(S) 1518 # Where v is velocity, n is manning's coefficient, h is depth and S is the slope into the domain. 1519 # Let mu be the momentum (vh), then this equation becomes: mu = 1/n h^{5/3} sqrt(S) 1520 # from which we can isolate depth to get 1521 # h = (mu n/sqrt(S) )^{3/5} 1522 1523 slope = 0 # get gradient for this triangle dot normal 1524 1525 # get manning coef from this triangle 1526 friction = self.domain.get_quantity('friction').get_values(location='edges', 1527 indices=[vol_id])[0] 1528 mannings_n = friction[edge_id] 1529 1530 if slope > epsilon and mannings_n > epsilon: 1531 depth = pow(self.average_momentum * mannings_n/math.sqrt(slope), 3.0/5) 1532 else: 1533 depth = 1.0 1534 1535 # Elevation on this edge 1536 1537 z = self.domain.get_quantity('elevation').get_values(location='edges', 1538 indices=[vol_id])[0] 1539 elevation = z[edge_id] 1540 1541 # Assign conserved quantities and return 1542 q = num.array([elevation + depth, xmomentum, ymomentum], num.Float) 1543 return q 1544 1545 1546 1547 1548 1549 1467 1550 class Field_boundary(Boundary): 1468 1551 """Set boundary from given field represented in an sww file containing
Note: See TracChangeset
for help on using the changeset viewer.