source: inundation/ga/storm_surge/examples/beach.py @ 668

"""Example of the inundationmodel.

A wave of water is washed up ontp a hypothetical beach.

To run:

python beach.py
"""

######################
# Module imports 
#

import sys
from os import sep, path
sys.path.append('..'+sep+'pyvolution')

from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
     Transmissive_boundary, Time_boundary, Wind_stress

from pmesh2domain import pmesh_to_domain_instance
from config import data_dir
from util import read_polygon, Polygon_function
from math import pi
from Numeric import choose, greater, ones, sin


######################
# Domain
name = 'beach'
print 'Creating domain from %s.tsh' %name
domain = pmesh_to_domain_instance(name + '.tsh', Domain)

#HACK to remove internal boundary (until we get a solution)
for b in domain.boundary.keys():
    if domain.boundary[b] == '':
        del domain.boundary[b]

domain = Domain(domain.coordinates, domain.triangles, domain.boundary)

domain.store = True
domain.set_name(name)
print "Output being written to " + data_dir + sep + \
      domain.filename + "_size%d." %len(domain) + domain.format


def bathymetry(x, y):
    #cut = 55
    #return -(x - cut)/5

    cut = 75
    return choose( greater(x, cut), ( -(x - 55)/5, -4*ones(x.shape) )) 
    ####return choose( greater(x, cut), (-20*ones(x.shape), ones(x.shape)))   
    

def topography(x, y):
    return 4*sin(pi*x/50) + 1

def riverbed(x, y):
    #return (y-100)/10 + 5*x/100
    return (y-100)/50 + 4.4
    

shoreline = [[40,0], [100,0], [100,100], [65,100], 
             [55,90], [55,70], [56,50], [50,25], [40,0]]        

land = [[65,100], [55,90], [55,70], [56,50], [50,25], [40,0],
        [0,0], [0,100]]        
             
building1 = [[45,80], [50,78], [52,83], [45,83]]             
building2 = [[35,75], [40,75], [40,80], [35,80]]             
building3 = [[42,63], [47,63], [50,68], [37,68]]             
building4 = [[28,70], [28,60], [33,60], [33,70]]             
building5 = [[10,70], [10,60], [15,60], [15,70]]             
building6 = [[10,60], [10,50], [15,50], [15,60]]             




river = [[20,100], [18,90], [20,80], [20,60], [15,40], [11,20], [2,0], [10,0],
         [14,10], [20,30], [24,45], [27,80], [27,85], [35,90], [39,100]]
              
print 'Set elevation'    
domain.set_quantity('elevation',
    Polygon_function( [(shoreline, bathymetry), (land, topography), 
                      (building1, 3), (building2, 9), 
                      (building3, 8), (building4, 7),
                      (building5, 8), (building6, 6),
                      (river, riverbed)]))
                      
print 'Set level'                     
domain.set_quantity('level', 
                    Polygon_function( [(shoreline, -1.5), 
                    (land, -10)] )) 
                    
print 'Set friction'                                        
domain.set_quantity('friction', 0.03)

print domain.get_extent()



#Add lateral wind gusts bearing 135 degrees
def gust(t,x,y): 
    from math import sin, pi
    from Numeric import zeros, ones, Float

    N = len(x)

    tt = sin(2*pi*t/50)

    if tt > 0.7:
        return 10000*tt*ones(N, Float)
    else:
        return zeros(N, Float)
    
domain.forcing_terms.append(Wind_stress(gust, 135))


######################
# Boundary conditions

print 'Boundaries'
Br = Reflective_boundary(domain)
Bo = Transmissive_boundary(domain)

#Constant inflow
Bd = Dirichlet_boundary([0.0, -1.0, 0.0])
Bt = Time_boundary(domain, lambda t: [ 3.0*(1+sin(2*pi*t/100)), 0.0, 0.0])

print 'Available boundary tags are', domain.get_boundary_tags()

#Set boundary conditions
tags = {}
tags['ocean'] = Bt
tags['wall'] = Br 
tags['outflow'] = Bo 
tags['exterior'] = Br
tags['external'] = Br 
domain.set_boundary(tags)


domain.check_integrity()

######################
#Evolution
for t in domain.evolve(yieldstep = 0.2, finaltime = 500):
    domain.write_time()
    
print 'Done'