Ignore:
Timestamp:
Dec 6, 2004, 1:02:40 AM (20 years ago)
Author:
ole
Message:

Showcase - beach example

File:
1 edited

Legend:

Unmodified
Added
Removed
  • inundation/ga/storm_surge/examples/beach.py

    r668 r672  
    2323from util import read_polygon, Polygon_function
    2424from math import pi
    25 from Numeric import choose, greater, ones, sin
    26 
     25from Numeric import choose, greater, ones, sin, exp
     26import time
    2727
    2828######################
     
    3333
    3434#HACK to remove internal boundary (until we get a solution)
     35#Is this still necessary?
    3536for b in domain.boundary.keys():
    3637    if domain.boundary[b] == '':
     
    4647
    4748def bathymetry(x, y):
    48     #cut = 55
    49     #return -(x - cut)/5
    50 
    5149    cut = 75
    52     return choose( greater(x, cut), ( -(x - 55)/5, -4*ones(x.shape) ))
    53     ####return choose( greater(x, cut), (-20*ones(x.shape), ones(x.shape)))   
     50   
     51    res = choose( greater(x, cut), ( -(x - 55)/5, -4*ones(x.shape) ))
     52    res += (100-y)/50 - 0.75
     53    return res
    5454   
    5555
    5656def topography(x, y):
    57     return 4*sin(pi*x/50) + 1
     57
     58    import RandomArray
     59    #z = 4*sin(pi*x/50) + (100-y)/50 + 1 #Ridge
     60   
     61    z = -4.0*x/25 + 8  + (100-y)/50 #Beach
     62   
     63    #z += sin(pi*x/5)/4
     64    #z += RandomArray.normal(0, 1, x.shape)/4
     65   
     66    z += 6*exp( -((x-30)/10)**2 ) * exp( -((y-50)/8)**2 )     #Mound
     67    z += 4*exp( -((x-30)/4)**2 ) * exp( -((y-26)/10)**2 )     #extra ridge   
     68    z += 4*exp( -((x-30)/5)**2 ) * exp( -((y-10)/8)**2 )     #extra ridge       
     69   
     70    return z
    5871
    5972def riverbed(x, y):
    60     #return (y-100)/10 + 5*x/100
    61     return (y-100)/50 + 4.4
     73    return (y-100)/70 - 4.0*x/25 + 8
    6274   
    6375
     
    6678
    6779land = [[65,100], [55,90], [55,70], [56,50], [50,25], [40,0],
    68         [0,0], [0,100]]       
     80        [0,0], [0,100]]
     81                       
    6982             
    70 building1 = [[45,80], [50,78], [52,83], [45,83]]             
     83building1 = [[45,80], [49,78], [52,83], [46,83]]             
    7184building2 = [[35,75], [40,75], [40,80], [35,80]]             
    72 building3 = [[42,63], [47,63], [50,68], [37,68]]             
    73 building4 = [[28,70], [28,60], [33,60], [33,70]]             
    74 building5 = [[10,70], [10,60], [15,60], [15,70]]             
    75 building6 = [[10,60], [10,50], [15,50], [15,60]]             
    76 
    77 
    78 
     85building3 = [[42,63], [46,61], [48,65], [44,67]]             
     86building4 = [[28,56], [28,51], [33,51], [33,56]]             
     87building5 = [[10,70], [10,65], [15,65], [15,70]]             
     88building6 = [[10,50], [10,45], [15,45], [15,50]]             
    7989
    8090river = [[20,100], [18,90], [20,80], [20,60], [15,40], [11,20], [2,0], [10,0],
     
    8292             
    8393print 'Set elevation'   
     94t0 = time.time()
    8495domain.set_quantity('elevation',
    8596    Polygon_function( [(shoreline, bathymetry), (land, topography),
    86                       (building1, 3), (building2, 9),
    87                       (building3, 8), (building4, 7),
    88                       (building5, 8), (building6, 6),
     97                      (building1, 6), (building2, 8),
     98                      (building3, 6), (building4, 13),
     99                      (building5, 10), (building6, 11),
    89100                      (river, riverbed)]))
     101                     
     102print 'That took %.2f seconds' %(time.time()-t0)
    90103                     
    91104print 'Set level'                     
     
    95108                   
    96109print 'Set friction'                                       
    97 domain.set_quantity('friction', 0.03)
     110domain.set_quantity('friction', 0.075)
    98111
    99112print domain.get_extent()
     
    111124
    112125    if tt > 0.7:
    113         return 10000*tt*ones(N, Float)
     126        return 30000*tt*ones(N, Float)
    114127    else:
    115128        return zeros(N, Float)
     
    125138Bo = Transmissive_boundary(domain)
    126139
    127 #Constant inflow
    128 Bd = Dirichlet_boundary([0.0, -1.0, 0.0])
    129 Bt = Time_boundary(domain, lambda t: [ 3.0*(1+sin(2*pi*t/100)), 0.0, 0.0])
     140#Constant outflow
     141Bd = Dirichlet_boundary([-10, 0.0, 0.0])
     142Bt = Time_boundary(domain, lambda t: [ 2.8*(1+sin(2*pi*t/100)), 0.0, 0.0])
    130143
    131144print 'Available boundary tags are', domain.get_boundary_tags()
     
    135148tags['ocean'] = Bt
    136149tags['wall'] = Br
    137 tags['outflow'] = Bo
     150tags['wall1'] = Br
     151tags['outflow'] = Bd
    138152tags['exterior'] = Br
    139 tags['external'] = Br
     153tags['external'] = Br
     154tags['land'] = Bo 
     155tags['westbank'] = None    #Riverbank
     156tags['eastbank'] = None    #Riverbank
     157tags['eastbankN'] = None    #Riverbank
    140158domain.set_boundary(tags)
    141159
     
    145163######################
    146164#Evolution
    147 for t in domain.evolve(yieldstep = 0.2, finaltime = 500):
     165for t in domain.evolve(yieldstep = 0.2, finaltime = 300):
    148166    domain.write_time()
    149167   
Note: See TracChangeset for help on using the changeset viewer.