source: documentation/user_manual/examples/harbour.py @ 3563

Last change on this file since 3563 was 3563, checked in by ole, 16 years ago

Moved shallow water out from the old pyvolution directory.
All tests pass, most examples and okushiri works again.

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