source: inundation/pyvolution/netherlands-chris.py @ 1882

Last change on this file since 1882 was 1290, checked in by steve, 19 years ago

visualisation stuff

File size: 6.7 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
99<<<<<<< .mine
100N = 4
101=======
102#N = 140
103>>>>>>> .r817
104
105#N = 15
106
107print 'Creating domain'
108#Create basic mesh
109points, vertices, boundary = rectangular(N, N)
110
111#Create shallow water domain
112domain = Domain(points, vertices, boundary)
113   
114domain.check_integrity()
115domain.default_order = 2
116#domain.beta_h=0
117
118#Output params
119domain.smooth = True
120domain.reduction = min  #Looks a lot better on top of steep slopes
121
122print "Number of triangles = ", len(domain)
123
124
125if N > 40:
126    domain.visualise = False
127    domain.checkpoint = False
128    domain.store = True    #Store for visualisation purposes
129    domain.format = 'sww'   #Native netcdf visualisation format
130    import sys, os
131    #FIXME: This was os.path.splitext but caused weird filenames based on root
132    base = os.path.basename(sys.argv[0])
133    domain.filename, _ = os.path.splitext(base)
134else:
135    domain.visualise = False
136    domain.checkpoint = False
137    domain.store = False   
138
139   
140#Set bed-slope and friction
141inflow_stage = 0.08
142manning = 0.02
143Z = Weir(inflow_stage)
144
145print 'Field values'
146domain.set_quantity('elevation', Z)
147domain.set_quantity('friction', manning)
148
149
150######################
151# Boundary conditions
152#
153print 'Boundaries'
154Br = Reflective_boundary(domain)
155Bt = Transmissive_boundary(domain)
156
157#Constant inflow
158Bd = Dirichlet_boundary([2*inflow_stage, 0.0, 0.0])
159
160
161#Set boundary conditions
162domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
163                 
164
165######################
166#Initial condition
167#
168print 'Initial condition'
169domain.set_quantity('stage', Constant_height(Z, 0.))
170
171#Evolve
172import time
173t0 = time.time()
174
175for t in domain.evolve(yieldstep = 0.5, finaltime = 1.0):
176    domain.write_time()
177
178
179
180   
181print 'That took %.2f seconds' %(time.time()-t0)
182print 'time', domain.write_time()
183
184print domain.coordinates
185print '*****'
186print domain.vertex_coordinates
187print '*****'
188print domain.quantities['xmomentum'].centroid_values
189print '*****'
190print domain.quantities['xmomentum'].edge_values
191print '*****'
192print domain.quantities['stage'].vertex_values
193print '*****'
194print domain.quantities['stage'].explicit_update
195
196from shallow_water import *
197
198def compute_fluxes_python(domain):
199    """Compute all fluxes and the timestep suitable for all volumes
200    in domain.
201
202    Compute total flux for each conserved quantity using "flux_function"
203
204    Fluxes across each edge are scaled by edgelengths and summed up
205    Resulting flux is then scaled by area and stored in
206    explicit_update for each of the three conserved quantities
207    stage, xmomentum and ymomentum   
208
209    The maximal allowable speed computed by the flux_function for each volume
210    is converted to a timestep that must not be exceeded. The minimum of
211    those is computed as the next overall timestep.
212
213    Post conditions:
214      domain.explicit_update is reset to computed flux values
215      domain.timestep is set to the largest step satisfying all volumes.
216    """
217
218    import sys
219    from Numeric import zeros, Float
220
221    N = domain.number_of_elements
222   
223    #Shortcuts
224    Stage = domain.quantities['stage']
225    Xmom = domain.quantities['xmomentum']
226    Ymom = domain.quantities['ymomentum']
227    Bed = domain.quantities['elevation']       
228
229    #Arrays
230    stage = Stage.edge_values
231    xmom =  Xmom.edge_values
232    ymom =  Ymom.edge_values
233    bed =   Bed.edge_values   
234
235    stage_bdry = Stage.boundary_values
236    xmom_bdry =  Xmom.boundary_values
237    ymom_bdry =  Ymom.boundary_values
238   
239    flux = zeros((N,3), Float) #Work array for summing up fluxes
240
241    #Loop
242    timestep = float(sys.maxint)   
243    for k in range(N):
244
245        for i in range(3):
246            #Quantities inside volume facing neighbour i
247            ql = [stage[k, i], xmom[k, i], ymom[k, i]]
248            zl = bed[k, i]
249
250            #Quantities at neighbour on nearest face
251            n = domain.neighbours[k,i] 
252            if n < 0:
253                m = -n-1 #Convert negative flag to index
254                qr = [stage_bdry[m], xmom_bdry[m], ymom_bdry[m]]
255                zr = zl #Extend bed elevation to boundary
256            else:   
257                m = domain.neighbour_edges[k,i]
258                qr = [stage[n, m], xmom[n, m], ymom[n, m]]
259                zr = bed[n, m]               
260
261               
262            #Outward pointing normal vector   
263            normal = domain.normals[k, 2*i:2*i+2]
264
265            #Flux computation using provided function
266            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
267       
268            flux[k,:] = edgeflux
269
270    return flux
271   
272flux = compute_fluxes_python(domain)
273print 'flux'
274print flux
275
276   
277# THis was pulled out of
Note: See TracBrowser for help on using the repository browser.