# source:anuga_core/source/anuga/examples/netherlands.py@3737

Last change on this file since 3737 was 3689, checked in by ole, 17 years ago

Added functionality for getting arbitrary interpolated values in Quantity as well as calculating inundation height and location. This work was done at SUT during the last week of September 2006.

File size: 5.2 KB
Line
1"""Example of shallow water wave equation.
2
3This is called Netherlands because it shows a dam with a gap in it and
4stylised housed behind it and below the water surface.
5
6"""
7
8######################
9# Module imports
10#
11#import rpdb
12#rpdb.set_active()
13
14from anuga.shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
15     Transmissive_boundary, Constant_height, Constant_stage
16
17from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
18from Numeric import array
19from anuga.visualiser.vtk_realtime_visualiser import Visualiser
20
21class Weir:
22    """Set a bathymetry for simple weir with a hole.
23    x,y are assumed to be in the unit square
24    """
25
26    def __init__(self, stage):
27        self.inflow_stage = stage
28
29    def __call__(self, x, y):
30        from Numeric import zeros, Float
31
32        N = len(x)
33        assert N == len(y)
34
35        z = zeros(N, Float)
36        for i in range(N):
37            z[i] = -x[i]/20  #General slope
38
39            #Flattish bit to the left
40            if x[i] <= 0.3:
41                #z[i] = -x[i]/5
42                z[i] = -x[i]/20
43
44
45            #Weir
46            if x[i] > 0.3 and x[i] < 0.4:
47                z[i] = -x[i]/20+1.2
48
49            #Dip
50            #if x[i] > 0.6 and x[i] < 0.9:
51            #    z[i] = -x[i]/20-0.5  #-y[i]/5
52
53            #Hole in weir
54            #if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
55            if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.4 and y[i] < 0.6:
56                #z[i] = -x[i]/5
57                z[i] = -x[i]/20
58
59            #Poles
60            #if x[i] > 0.65 and x[i] < 0.8 and y[i] > 0.55 and y[i] < 0.65 or\
61            #       x[i] > 0.75 and x[i] < 0.9 and y[i] > 0.35 and y[i] < 0.45:
62            #    z[i] = -x[i]/20+0.4
63
64            if (x[i] - 0.72)**2 + (y[i] - 0.6)**2 < 0.05**2:# or\
65                   #x[i] > 0.75 and x[i] < 0.9 and y[i] > 0.35 and y[i] < 0.45:
66                z[i] = -x[i]/20+0.4
67
68
69            #Wall
70            if x[i] > 0.995:
71                z[i] = -x[i]/20+0.3
72
73        return z/2
74
75
76
77######################
78# Domain
79#
80
81
82N = 150 #size = 45000
83N = 130 #size = 33800
84N = 600 #Size = 720000
85N = 100
86
87
88#N = 15
89
90print 'Creating domain'
91#Create basic mesh
92points, elements, boundary = rectangular_cross(N, N)
93
94#Create shallow water domain
95domain = Domain(points, elements, boundary, use_inscribed_circle=True)
96
97domain.check_integrity()
98
99#Setup order and all the beta's for the limiters (these should become defaults
100domain.default_order = 2
101domain.beta_w      = 1.0
102domain.beta_w_dry  = 0.2
103domain.beta_uh     = 1.0
104domain.beta_uh_dry = 0.2
105domain.beta_vh     = 1.0
106domain.beta_vh_dry = 0.2
107
108#Output params
109domain.smooth = False
110domain.reduction = min  #Looks a lot better on top of steep slopes
111
112print "Number of triangles = ", len(domain)
113print "Extent = ", domain.get_extent()
114
115#Set bed-slope and friction
116inflow_stage = 0.5
117manning = 0.03
118Z = Weir(inflow_stage)
119
120if N > 150:
121    domain.visualise = False
122    domain.checkpoint = False
123    domain.store = True    #Store for visualisation purposes
124    domain.format = 'sww'   #Native netcdf visualisation format
125    import sys, os
126    #FIXME: This was os.path.splitext but caused weird filenames based on root
127    base = os.path.basename(sys.argv[0])
128    domain.filename, _ = os.path.splitext(base)
129else:
130    #domain.initialise_visualiser(rect=[0.0,0.0,1.0,1.0])
131    #domain.initialise_visualiser()
132    #domain.visualiser.coloring['stage'] = False
133    domain.visualise_timer = True
134    domain.checkpoint = False
135    domain.store = False
136    vis = Visualiser(domain,title="netherlands")
137    vis.setup['elevation'] = True
138    vis.updating['stage'] = True
139    vis.qcolor['stage'] = (0.0,0.0,0.8)
140    vis.coloring['stage']= True
141
142
143
144print 'Field values'
145domain.set_quantity('elevation', Z)
146domain.set_quantity('friction', manning)
147
148
149######################
150# Boundary conditions
151#
152print 'Boundaries'
153Br = Reflective_boundary(domain)
154Bt = Transmissive_boundary(domain)
155
156#Constant inflow
157Bd = Dirichlet_boundary([inflow_stage, 0.0, 0.0])
158
159
160#Set boundary conditions
161domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
162
163
164######################
165#Initial condition
166#
167print 'Initial condition'
168domain.set_quantity('stage', expression='Z + 0.0')
169
170#Evolve
171import time
172t0 = time.time()
173
174
175#from realtime_visualisation_new import Visualiser
176#V = Visualiser(domain,title='netherlands')
177#V.update_quantity('stage')
178#V.update_quantity('elevation')
179
180
181
182
183for t in domain.evolve(yieldstep = 0.005, finaltime = None):
184    domain.write_time()
185    #domain.write_boundary()
186    vis.update()
187    print domain.quantities['stage'].get_values(location='centroids',
188                                                indices=[0])
189    #domain.visualiser.update_quantity('elevation')
190    #time.sleep(0.1)
191    #raw_input('pause>')
192    #V.update_quantity('stage')
193    #rpdb.set_active()
194    #domain.visualiser.scale_z = 1.0
195    #domain.visualiser.update_quantity_color('xmomentum',scale_z = 4.0)
196    #integral_label.text='Integral=%10.5e'%domain.quantities['stage'].get_integral()
197
198
199print 'That took %.2f seconds' %(time.time()-t0)
200vis.shutdown()
Note: See TracBrowser for help on using the repository browser.