Changeset 7348


Ignore:
Timestamp:
Aug 10, 2009, 5:40:57 AM (15 years ago)
Author:
ole
Message:

Elaborated on variable bedslope example for use with AOGS2009 talk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/documentation/user_manual/demos/channel_variable.py

    r7346 r7348  
    1616# Setup computational domain
    1717#------------------------------------------------------------------------------
    18 length = 40.
     18length = 24.
    1919width = 5.
    20 dx = dy = 1 #.1           # Resolution: Length of subdivisions on both axes
     20dx = dy = 0.05 #.1           # Resolution: Length of subdivisions on both axes
    2121
    2222points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
    2323                                               len1=length, len2=width)
    2424domain = Domain(points, vertices, boundary)
    25 domain.set_name('channel_variable')                  # Output name
     25domain.set_name('channel_variable_2')                  # Output name
    2626print domain.statistics()
    2727domain.set_quantities_to_be_stored({'elevation': 1, #2,
     
    3535
    3636    z = -x/100
    37 
     37   
    3838    N = len(x)
    3939    for i in range(N):
    4040        # Step
    41         if 10 < x[i] < 12:
     41        if 2 < x[i] < 4:
    4242            z[i] += 0.4 - 0.05*y[i]
    43 
    44         # Constriction
    45         #if 27 < x[i] < 29 and y[i] > 3:
    46         #    z[i] += 2
    47 
    48         # Pole
    49         #if (x[i] - 34)**2 + (y[i] - 2)**2 < 0.4**2:
    50         #    z[i] += 2
    51 
     43   
     44        # Permanent pole
     45        if (x[i] - 8)**2 + (y[i] - 2)**2 < 0.4**2:
     46            z[i] += 1
    5247    return z
     48   
    5349
    5450def pole_increment(x,y):
     
    6157    N = len(x)
    6258    for i in range(N):
    63         # Pole
    64         if (x[i] - 20)**2 + (y[i] - 2)**2 < 0.4**2:
    65             z[i] += 0.05
     59        # Pole 1
     60        if (x[i] - 12)**2 + (y[i] - 3)**2 < 0.4**2:
     61            z[i] += 0.01
     62           
     63    for i in range(N):
     64        # Pole 2
     65        if (x[i] - 14)**2 + (y[i] - 2)**2 < 0.4**2:
     66            z[i] += 0.005
     67                       
    6668    return z
    6769   
     
    8789shrinking = False
    8890done = False
    89 for t in domain.evolve(yieldstep=0.1, finaltime=30.0):
     91for t in domain.evolve(yieldstep=0.1, finaltime=40.0):
    9092    print domain.timestepping_statistics()
    9193
     
    9496    #print 'Level at gauge point = %.2fm' % w
    9597           
    96     z = domain.get_quantity('elevation').\
    97         get_values(interpolation_points=[[20, 2]])           
    98     print 'Elevation at pole location = %.2fm' % z           
     98    #z = domain.get_quantity('elevation').\
     99    #    get_values(interpolation_points=[[12, 3]])           
     100    #print 'Elevation at pole location = %.2fm' % z           
    99101
    100102    # Start variable elevation after 10 seconds   
     
    102104        growing = True
    103105           
    104     # Turn around when pole has reached a height of 2 m
    105     if z >= 2 and growing:
     106    # Stop growing when pole has reached a certain height
     107    if t > 16 and growing:
    106108        growing = False
     109        shrinking = False
     110       
     111    # Start shrinking
     112    if t > 20:
    107113        shrinking = True
     114        growing = False
    108115       
    109     # Stop changing when pole has shrunk to 1 m
    110     if z <= 1 and shrinking:
     116    # Stop changing when pole has shrunk to original level
     117    if t > 25 and shrinking:
    111118        done = True
    112119        shrinking = growing = False
     120        domain.set_quantity('elevation', topography)
    113121
    114122    # Grow or shrink               
     
    117125       
    118126    if shrinking:   
    119         domain.add_quantity('elevation', lambda x,y: -pole_increment(x,y))   
     127        domain.add_quantity('elevation', lambda x,y: -2*pole_increment(x,y))   
    120128       
Note: See TracChangeset for help on using the changeset viewer.