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

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

Moved obsolete code and examples to their appropriate locations

File size: 6.6 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 = 4
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#Evolve
168import time
169t0 = time.time()
170
171for t in domain.evolve(yieldstep = 0.5, finaltime = 1.0):
172    domain.write_time()
173
174
175
176   
177print 'That took %.2f seconds' %(time.time()-t0)
178print 'time', domain.write_time()
179
180print domain.coordinates
181print '*****'
182print domain.vertex_coordinates
183print '*****'
184print domain.quantities['xmomentum'].centroid_values
185print '*****'
186print domain.quantities['xmomentum'].edge_values
187print '*****'
188print domain.quantities['stage'].vertex_values
189print '*****'
190print domain.quantities['stage'].explicit_update
191
192from shallow_water import *
193
194def compute_fluxes_python(domain):
195    """Compute all fluxes and the timestep suitable for all volumes
196    in domain.
197
198    Compute total flux for each conserved quantity using "flux_function"
199
200    Fluxes across each edge are scaled by edgelengths and summed up
201    Resulting flux is then scaled by area and stored in
202    explicit_update for each of the three conserved quantities
203    stage, xmomentum and ymomentum   
204
205    The maximal allowable speed computed by the flux_function for each volume
206    is converted to a timestep that must not be exceeded. The minimum of
207    those is computed as the next overall timestep.
208
209    Post conditions:
210      domain.explicit_update is reset to computed flux values
211      domain.timestep is set to the largest step satisfying all volumes.
212    """
213
214    import sys
215    from Numeric import zeros, Float
216
217    N = domain.number_of_elements
218   
219    #Shortcuts
220    Stage = domain.quantities['stage']
221    Xmom = domain.quantities['xmomentum']
222    Ymom = domain.quantities['ymomentum']
223    Bed = domain.quantities['elevation']       
224
225    #Arrays
226    stage = Stage.edge_values
227    xmom =  Xmom.edge_values
228    ymom =  Ymom.edge_values
229    bed =   Bed.edge_values   
230
231    stage_bdry = Stage.boundary_values
232    xmom_bdry =  Xmom.boundary_values
233    ymom_bdry =  Ymom.boundary_values
234   
235    flux = zeros((N,3), Float) #Work array for summing up fluxes
236
237    #Loop
238    timestep = float(sys.maxint)   
239    for k in range(N):
240
241        for i in range(3):
242            #Quantities inside volume facing neighbour i
243            ql = [stage[k, i], xmom[k, i], ymom[k, i]]
244            zl = bed[k, i]
245
246            #Quantities at neighbour on nearest face
247            n = domain.neighbours[k,i] 
248            if n < 0:
249                m = -n-1 #Convert negative flag to index
250                qr = [stage_bdry[m], xmom_bdry[m], ymom_bdry[m]]
251                zr = zl #Extend bed elevation to boundary
252            else:   
253                m = domain.neighbour_edges[k,i]
254                qr = [stage[n, m], xmom[n, m], ymom[n, m]]
255                zr = bed[n, m]               
256
257               
258            #Outward pointing normal vector   
259            normal = domain.normals[k, 2*i:2*i+2]
260
261            #Flux computation using provided function
262            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
263       
264            flux[k,:] = edgeflux
265
266    return flux
267   
268flux = compute_fluxes_python(domain)
269print 'flux'
270print flux
271
272   
273# THis was pulled out of
Note: See TracBrowser for help on using the repository browser.