Changeset 6644


Ignore:
Timestamp:
Mar 27, 2009, 2:54:08 PM (15 years ago)
Author:
ted
Message:

Added function to compute depth dependent friction.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r6636 r6644  
    14861486                xmom_update[k] += S*uh[k]
    14871487                ymom_update[k] += S*vh[k]
     1488
     1489
     1490
     1491
     1492def depth_dependent_friction(domain, default_friction,
     1493                             surface_roughness_data,
     1494                             verbose=False):
     1495    """Returns an array of friction values for each wet element adjusted for depth.
     1496
     1497    Inputs:
     1498        domain - computational domain object
     1499        default_friction - depth independent bottom friction
     1500        surface_roughness_data - N x 5 array of n0, d1, n1, d2, n2 values for each
     1501        friction region.
     1502
     1503    Outputs:
     1504        wet_friction - Array that can be used directly to update friction as follows:
     1505                       domain.set_quantity('friction', wet_friction)
     1506
     1507       
     1508       
     1509    """
     1510
     1511    import Numeric as num
     1512   
     1513    # Create a temp array to store updated depth dependent friction for wet elements
     1514    # EHR this is outwardly inneficient but not obvious how to avoid recreating each call??????
     1515    N=len(domain)
     1516    wet_friction    = num.zeros(N, num.Float)
     1517    wet_friction[:] = default_n0   # Initially assign default_n0 to all array so sure have no zeros values
     1518   
     1519   
     1520    depth = domain.create_quantity_from_expression('stage - elevation')  # create depth instance for this timestep
     1521    # Recompute depth as vector 
     1522    d = depth.get_values(location='centroids')
     1523 
     1524    # rebuild the 'friction' values adjusted for depth at this instant
     1525    for i in domain.get_wet_elements():                                  # loop for each wet element in domain
     1526       
     1527        # Get roughness data for each element
     1528        n0 = float(surface_roughness_data[i,0])
     1529        d1 = float(surface_roughness_data[i,1])
     1530        n1 = float(surface_roughness_data[i,2])
     1531        d2 = float(surface_roughness_data[i,3])
     1532        n2 = float(surface_roughness_data[i,4])
     1533       
     1534       
     1535        # Recompute friction values from depth for this element
     1536               
     1537        if d[i]   <= d1:
     1538            depth_dependent_friction = n1
     1539        elif d[i] >= d2:
     1540            depth_dependent_friction = n2
     1541        else:
     1542            depth_dependent_friction = n1+((n2-n1)/(d2-d1))*(d[i]-d1)
     1543           
     1544        # check sanity of result
     1545        if (depth_dependent_friction  < 0.010 or depth_dependent_friction > 9999.0) :
     1546            print model_data.basename+' >>>> WARNING: computed depth_dependent friction out of range ddf,n1,n2 ', depth_dependent_friction, n1,n2
     1547       
     1548        # update depth dependent friction  for that wet element
     1549        wet_friction[i] = depth_dependent_friction
     1550       
     1551    # EHR add code to show range of 'friction across domain at this instant as sanity check?????????
     1552   
     1553    if verbose :
     1554        nvals=domain.get_quantity('friction').get_values(location='centroids')        # return array of domain nvals
     1555        n_min=min(nvals)
     1556        n_max=max(nvals)
     1557       
     1558        print "         ++++ calculate_depth_dependent_friction  - Updated friction - range  %7.3f to %7.3f" %(n_min,n_max)
     1559   
     1560    return wet_friction
     1561
     1562
     1563
     1564
     1565
     1566
     1567
     1568
     1569
    14881570
    14891571
Note: See TracChangeset for help on using the changeset viewer.