source: trunk/anuga_core/documentation/user_manual/examples/harbour.py @ 7808

Last change on this file since 7808 was 7808, checked in by hudson, 12 years ago

Examples work with new API.

File size: 5.5 KB
Line 
1"""Simple water flow example using ANUGA
2
3Water driven towards a continental slope and into a harbour
4environment. Set up as per:
5P.L.-F Liu
6Effects of the Continental Shelf on Harbor Resonance, in
7Tsunamis - Their Science and Engineering,
8edited by K. Iida and T. Iwasaki, 303-314.
91983, Terra Scientific Publishing Company, Tokyo.
10
11"""
12
13#------------------------------------------------------------------------------
14# Import necessary modules
15#------------------------------------------------------------------------------
16import anuga
17
18waveheight = 1.0
19harbour_length = 2000.0 
20shelf_length = 10000.0
21bottom_length = 10000.0
22harbour_width = 250.0
23shelf_depth = -100. 
24harbour_depth = -10.
25if harbour_depth > shelf_depth: print 'WARNING: depths not consistent'
26#------------------------------------------------------------------------------
27# Setup computational domain
28#------------------------------------------------------------------------------
29
30east = 10000.+harbour_length+shelf_length+bottom_length
31west = 0.
32south = 0.
33north = 5000.
34polygon = [[east,north],[west,north],[west,south],[east,south]]
35meshname = 'harbour_test.msh'
36
37anuga.create_mesh_from_regions(polygon,
38                         boundary_tags={'top': [0],
39                                        'left': [1],
40                                        'bottom': [2],
41                                        'right': [3]},
42                         maximum_triangle_area=50000,
43                         filename=meshname)
44                         #interior_regions=interior_regions)
45
46domain = anuga.Domain(meshname,use_cache=True,verbose = True)
47domain.set_name('harbour')                  # Output to harbour.sww
48domain.set_datadir('.')                     # Use current directory for output
49domain.set_quantities_to_be_stored(['stage',# Store all conserved quantities
50                                    'xmomentum',
51                                    'ymomentum'])   
52
53
54#------------------------------------------------------------------------------
55# Setup initial conditions
56#------------------------------------------------------------------------------
57# if can't do this, then can set up x = arange(west,east,xstep)
58# and y = arange(south,north,ystep) and use Geospatial data to set elevation
59
60def topography(x,y):
61    from numpy import zeros, float
62    z = zeros(len(x), float)
63   
64    xcentre = (west+east)/2.
65    ycentre = (south+north)/2.
66    ycentreup = ycentre+harbour_width/2.
67    ycentredown = ycentre-harbour_width/2.
68    for i, xi in enumerate(x):
69        yi = y[i]
70        if xi >= xcentre: z[i] = shelf_depth
71        if xi > west and xi <= xcentre-harbour_length: 
72            if yi < ycentredown: z[i] = 0.0
73            if yi > ycentreup: z[i] = 0.0
74            if yi > ycentredown and yi < ycentreup: z[i] = harbour_depth
75        if xi > xcentre-harbour_length and xi < xcentre: z[i] = harbour_depth
76       
77    return z
78
79domain.set_quantity('elevation', topography) # Use function for elevation
80domain.set_quantity('friction', 0. )         # Constant friction
81domain.set_quantity('stage', 0.0)            # Constant initial stage
82
83#------------------------------------------------------------------------------
84# Setup boundary conditions
85#------------------------------------------------------------------------------
86
87from math import sin, pi
88Br = anuga.Reflective_boundary(domain)      # Solid reflective wall
89Bt = anuga.Transmissive_boundary(domain)    # Continue all values on boundary
90Bd = anuga.Dirichlet_boundary([0.,0.,0.])   # Constant boundary values
91Bw = anuga.Time_boundary(domain=domain,     # Time dependent boundary 
92                   f=lambda t: [((10<t<20)*waveheight), 0.0, 0.0])
93
94# Associate boundary tags with boundary objects
95domain.set_boundary({'left': Br, 'right': Bw, 'top': Br, 'bottom': Br})
96
97#------------------------------------------------------------------------------
98# Evolve system through time
99#------------------------------------------------------------------------------
100
101stagestep = []
102for t in domain.evolve(yieldstep = 10, finaltime = 100.0):
103    domain.write_time()
104
105#-----------------------------------------------------------------------------
106# Interrogate solution
107#-----------------------------------------------------------------------------
108
109# Gauge line
110def gauge_line(west,east,north,south):
111    from numpy import arange
112    gaugex = arange(east,west,-1000.)
113    gauges = []
114    gaugey = []
115    for x in gaugex:
116        y = (north+south)/2.
117        gaugey.append(y)
118        gauges.append([x,y])
119    return gauges, gaugex, gaugey
120
121gauges, gaugex, gaugey = gauge_line(west,east,north,south)
122
123f = anuga.file_function('harbour.sww',
124                  quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum'],
125                  interpolation_points = gauges,
126                  verbose = True,
127                  use_cache = True)
128
129from pylab import plot, xlabel, ylabel, title, ion, close, savefig, figure
130ion()
131
132maxw = []
133minw = []
134maxd = []
135count = 0
136for k, g in enumerate(gauges):
137    stage = []
138    elevation = []
139    depths = []
140    for i, t in enumerate(f.get_time()):
141        w = f(t, point_id = k)[0]
142        elev = f(t, point_id = k)[1]
143        depth = w-elev
144        stage.append(w)
145        elevation.append(elev)
146        depths.append(elev)
147    maxw.append(max(stage))
148    minw.append(min(stage))
149    maxd.append(max(depths))
150
151filename = 'point'+str(waveheight)+'.csv'
152fid = open(filename,'w')   
153s = 'Waveheight,\n'
154fid.write(s)
155
156figure(1)
157plot(gaugex[:loc],maxw[:loc],'g-')
158xlabel('x')
159ylabel('stage')
160title('Maximum stage for gauge line')
161savefig('max_stage')
162
163close('all')
Note: See TracBrowser for help on using the repository browser.