source: anuga_core/source/anuga/examples/netherlands.py @ 3678

Last change on this file since 3678 was 3678, checked in by steve, 17 years ago

Added more limiting to cells near dry cells, use beta_*_dry

File size: 5.1 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', Constant_height(Z, 0.0))
169#domain.set_quantity('stage', Constant_stage(inflow_stage/2.0))
170
171#Evolve
172import time
173t0 = time.time()
174
175
176#from realtime_visualisation_new import Visualiser
177#V = Visualiser(domain,title='netherlands')
178#V.update_quantity('stage')
179#V.update_quantity('elevation')
180
181
182
183
184for t in domain.evolve(yieldstep = 0.005, finaltime = None):
185    domain.write_time()
186    #domain.write_boundary()
187    vis.update()
188    print domain.quantities['stage'].get_values(location='centroids',
189                                                indices=[0])
190    #domain.visualiser.update_quantity('elevation')
191    #time.sleep(0.1)
192    #raw_input('pause>')
193    #V.update_quantity('stage')
194    #rpdb.set_active()
195    #domain.visualiser.scale_z = 1.0
196    #domain.visualiser.update_quantity_color('xmomentum',scale_z = 4.0)
197    #integral_label.text='Integral=%10.5e'%domain.quantities['stage'].get_integral()
198
199
200print 'That took %.2f seconds' %(time.time()-t0)
201vis.shutdown()
Note: See TracBrowser for help on using the repository browser.