source: anuga_core/source/obsolete_code/netherlands-chris.py @ 5847

Last change on this file since 5847 was 5847, checked in by steve, 15 years ago

Changed parallel_api so that global mesh only needs to
be constructed on processor 0

File size: 7.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#
11from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
12     Transmissive_boundary, Constant_height
13
14from mesh_factory import rectangular
15from Numeric import array
16   
17
18class Weir:
19    """Set a bathymetry for simple weir with a hole.
20    x,y are assumed to be in the unit square
21    """
22   
23    def __init__(self, stage):
24        self.inflow_stage = stage
25
26    def __call__(self, x, y):   
27        from Numeric import zeros, Float
28   
29        N = len(x)
30        assert N == len(y)   
31
32        z = zeros(N, Float)
33        for i in range(N):
34            z[i] = -x[i]/20  #General slope
35
36            #Flattish bit to the left
37            if x[i] <= 0.3:
38                #z[i] = -x[i]/5
39                z[i] = -x[i]/20
40               
41           
42            #Weir
43            if x[i] > 0.3 and x[i] < 0.4:
44                z[i] = -x[i]/20+1.2
45
46            #Dip
47            #if x[i] > 0.6 and x[i] < 0.9:
48            #    z[i] = -x[i]/20-0.5  #-y[i]/5
49
50            #Hole in weir
51            #if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
52            if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.4 and y[i] < 0.6:
53                #z[i] = -x[i]/5
54                z[i] = -x[i]/20
55           
56            #Poles
57            #if x[i] > 0.65 and x[i] < 0.8 and y[i] > 0.55 and y[i] < 0.65 or\
58            #       x[i] > 0.75 and x[i] < 0.9 and y[i] > 0.35 and y[i] < 0.45:
59            #    z[i] = -x[i]/20+0.4
60
61            if (x[i] - 0.72)**2 + (y[i] - 0.6)**2 < 0.05**2:# or\
62                   #x[i] > 0.75 and x[i] < 0.9 and y[i] > 0.35 and y[i] < 0.45:
63                z[i] = -x[i]/20+0.4
64
65            #Wall
66            if x[i] > 0.995:
67                z[i] = -x[i]/20+0.3               
68
69        return z/2
70       
71   
72   
73######################
74# Domain
75#
76
77N = 250
78#N= 8
79N = 16
80#N = 4
81#N = 102
82N = 25
83N = 16
84N = 60
85N = 150 #size = 45000
86N = 130 #size = 33800
87#N = 60
88#N = 40
89N = 260
90#N = 150
91N = 264
92
93N = 600 #Size = 720000
94N = 20
95#N = 150
96N = 110
97N = 60
98
99N = 40
100#N = 140
101#N = 15
102
103print 'Creating domain'
104#Create basic mesh
105points, vertices, boundary = rectangular(N, N)
106
107#Create shallow water domain
108domain = Domain(points, vertices, boundary)
109   
110domain.check_integrity()
111domain.default_order = 2
112#domain.beta_h=0
113
114#Output params
115domain.smooth = True
116domain.reduction = min  #Looks a lot better on top of steep slopes
117
118print "Number of triangles = ", len(domain)
119
120
121if N > 40:
122    domain.visualise = False
123    domain.checkpoint = False
124    domain.store = True    #Store for visualisation purposes
125    domain.format = 'sww'   #Native netcdf visualisation format
126    import sys, os
127    #FIXME: This was os.path.splitext but caused weird filenames based on root
128    base = os.path.basename(sys.argv[0])
129    domain.filename, _ = os.path.splitext(base)
130else:
131    domain.visualise = False
132    domain.checkpoint = False
133    domain.store = False   
134
135   
136#Set bed-slope and friction
137inflow_stage = 0.08
138manning = 0.02
139Z = Weir(inflow_stage)
140
141print 'Field values'
142domain.set_quantity('elevation', Z)
143domain.set_quantity('friction', manning)
144
145
146######################
147# Boundary conditions
148#
149print 'Boundaries'
150Br = Reflective_boundary(domain)
151Bt = Transmissive_boundary(domain)
152
153#Constant inflow
154Bd = Dirichlet_boundary([2*inflow_stage, 0.0, 0.0])
155
156
157#Set boundary conditions
158domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
159                 
160
161######################
162#Initial condition
163#
164print 'Initial condition'
165domain.set_quantity('stage', Constant_height(Z, 0.))
166
167
168visualize = True
169if visualize:
170    from anuga.visualiser import RealtimeVisualiser
171    vis = RealtimeVisualiser(domain)
172    vis.render_quantity_height("elevation", zScale=100, offset = 5.0, dynamic=False)
173    vis.render_quantity_height("stage", zScale =100, dynamic=True, opacity = 0.3)
174    vis.colour_height_quantity('stage', (lambda q:q['stage'], -0.5, 0.5))
175    vis.start()
176
177time.sleep(2.0)
178
179
180
181#Evolve
182import time
183t0 = time.time()
184
185for t in domain.evolve(yieldstep = 0.5, finaltime = 1.0):
186    domain.write_time()
187
188    if visualize: vis.update()
189   
190if visualize: vis.evolveFinished()
191
192   
193print 'That took %.2f seconds' %(time.time()-t0)
194print 'time', domain.write_time()
195
196print domain.coordinates
197print '*****'
198print domain.vertex_coordinates
199print '*****'
200print domain.quantities['xmomentum'].centroid_values
201print '*****'
202print domain.quantities['xmomentum'].edge_values
203print '*****'
204print domain.quantities['stage'].vertex_values
205print '*****'
206print domain.quantities['stage'].explicit_update
207
208from shallow_water import *
209
210def compute_fluxes_python(domain):
211    """Compute all fluxes and the timestep suitable for all volumes
212    in domain.
213
214    Compute total flux for each conserved quantity using "flux_function"
215
216    Fluxes across each edge are scaled by edgelengths and summed up
217    Resulting flux is then scaled by area and stored in
218    explicit_update for each of the three conserved quantities
219    stage, xmomentum and ymomentum   
220
221    The maximal allowable speed computed by the flux_function for each volume
222    is converted to a timestep that must not be exceeded. The minimum of
223    those is computed as the next overall timestep.
224
225    Post conditions:
226      domain.explicit_update is reset to computed flux values
227      domain.timestep is set to the largest step satisfying all volumes.
228    """
229
230    import sys
231    from Numeric import zeros, Float
232
233    N = domain.number_of_elements
234   
235    #Shortcuts
236    Stage = domain.quantities['stage']
237    Xmom = domain.quantities['xmomentum']
238    Ymom = domain.quantities['ymomentum']
239    Bed = domain.quantities['elevation']       
240
241    #Arrays
242    stage = Stage.edge_values
243    xmom =  Xmom.edge_values
244    ymom =  Ymom.edge_values
245    bed =   Bed.edge_values   
246
247    stage_bdry = Stage.boundary_values
248    xmom_bdry =  Xmom.boundary_values
249    ymom_bdry =  Ymom.boundary_values
250   
251    flux = zeros((N,3), Float) #Work array for summing up fluxes
252
253    #Loop
254    timestep = float(sys.maxint)   
255    for k in range(N):
256
257        for i in range(3):
258            #Quantities inside volume facing neighbour i
259            ql = [stage[k, i], xmom[k, i], ymom[k, i]]
260            zl = bed[k, i]
261
262            #Quantities at neighbour on nearest face
263            n = domain.neighbours[k,i] 
264            if n < 0:
265                m = -n-1 #Convert negative flag to index
266                qr = [stage_bdry[m], xmom_bdry[m], ymom_bdry[m]]
267                zr = zl #Extend bed elevation to boundary
268            else:   
269                m = domain.neighbour_edges[k,i]
270                qr = [stage[n, m], xmom[n, m], ymom[n, m]]
271                zr = bed[n, m]               
272
273               
274            #Outward pointing normal vector   
275            normal = domain.normals[k, 2*i:2*i+2]
276
277            #Flux computation using provided function
278            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
279       
280            flux[k,:] = edgeflux
281
282    return flux
283   
284flux = compute_fluxes_python(domain)
285print 'flux'
286print flux
287
288   
289# THis was pulled out of
Note: See TracBrowser for help on using the repository browser.