# source:inundation/ga/storm_surge/pyvolution/netherlands.py@1290

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

visualisation stuff

File size: 3.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
77
78N = 150 #size = 45000
79N = 130 #size = 33800
80N = 600 #Size = 720000
81N = 30
82
83#N = 15
84
85print 'Creating domain'
86#Create basic mesh
87points, vertices, boundary = rectangular(N, N)
88
89#Create shallow water domain
90domain = Domain(points, vertices, boundary)
91
92domain.check_integrity()
93domain.default_order = 2
94#domain.beta_h=0
95
96#Output params
97domain.smooth = False
98domain.reduction = min  #Looks a lot better on top of steep slopes
99
100print "Number of triangles = ", len(domain)
101
102
103if N > 150:
104    domain.visualise = False
105    domain.checkpoint = False
106    domain.store = True    #Store for visualisation purposes
107    domain.format = 'sww'   #Native netcdf visualisation format
108    import sys, os
109    #FIXME: This was os.path.splitext but caused weird filenames based on root
110    base = os.path.basename(sys.argv[0])
111    domain.filename, _ = os.path.splitext(base)
112else:
113    domain.visualise = True
114    domain.checkpoint = False
115    domain.store = False
116
117
118#Set bed-slope and friction
119inflow_stage = 0.1
120manning = 0.02
121Z = Weir(inflow_stage)
122
123print 'Field values'
124domain.set_quantity('elevation', Z)
125domain.set_quantity('friction', manning)
126
127
128######################
129# Boundary conditions
130#
131print 'Boundaries'
132Br = Reflective_boundary(domain)
133Bt = Transmissive_boundary(domain)
134
135#Constant inflow
136Bd = Dirichlet_boundary([2*inflow_stage, 0.0, 0.0])
137
138
139#Set boundary conditions
140domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
141
142
143######################
144#Initial condition
145#
146print 'Initial condition'
147domain.set_quantity('stage', Constant_height(Z, 0.))
148
149#Evolve
150import time
151t0 = time.time()
152
153for t in domain.evolve(yieldstep = 0.2, finaltime = 10.0):
154    domain.write_time()
155
156
157print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.