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

Last change on this file since 3483 was 3483, checked in by ole, 18 years ago

Fixed topography for harbour example.

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 pmesh.mesh_interface import create_mesh_from_regions
18from pyvolution.shallow_water import Domain
19from pyvolution.shallow_water import Reflective_boundary
20from pyvolution.shallow_water import Dirichlet_boundary
21from pyvolution.shallow_water import Time_boundary
22from pyvolution.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 pyvolution.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.