1 | #! /usr/bin/python |
---|
2 | |
---|
3 | |
---|
4 | """ |
---|
5 | Forcing terms for the shallow water equations, gravity, friction etc |
---|
6 | """ |
---|
7 | |
---|
8 | __author__="Stephen Roberts" |
---|
9 | __date__ ="$05/06/2010 5:49:35 PM$" |
---|
10 | |
---|
11 | |
---|
12 | def gravity_for_loops(domain): |
---|
13 | """Apply gravitational pull in the presence of bed slope |
---|
14 | """ |
---|
15 | |
---|
16 | from anuga_1d.base.util import gradient |
---|
17 | |
---|
18 | xmom = domain.quantities['xmomentum'].explicit_update |
---|
19 | |
---|
20 | Stage = domain.quantities['stage'] |
---|
21 | Elevation = domain.quantities['elevation'] |
---|
22 | |
---|
23 | h = Stage.vertex_values - Elevation.vertex_values |
---|
24 | b = Elevation.vertex_values |
---|
25 | w = Stage.vertex_values |
---|
26 | |
---|
27 | x = domain.get_vertex_coordinates() |
---|
28 | g = domain.g |
---|
29 | |
---|
30 | for k in range(domain.number_of_elements): |
---|
31 | avg_h = sum( h[k,:] )/2 |
---|
32 | |
---|
33 | #Compute bed slope |
---|
34 | #x0, y0, x1, y1, x2, y2 = x[k,:] |
---|
35 | x0, x1 = x[k,:] |
---|
36 | #z0, z1, z2 = v[k,:] |
---|
37 | b0, b1 = b[k,:] |
---|
38 | |
---|
39 | |
---|
40 | #zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2) |
---|
41 | bx = gradient(x0, x1, b0, b1) |
---|
42 | |
---|
43 | #Update momentum (explicit update is reset to source values) |
---|
44 | xmom[k] += -g*bx*avg_h |
---|
45 | #xmom[k] = -g*bx*avg_h |
---|
46 | #stage[k] = 0.0 |
---|
47 | |
---|
48 | |
---|
49 | def gravity(domain): |
---|
50 | """Apply gravitational pull in the presence of bed slope |
---|
51 | """ |
---|
52 | |
---|
53 | xmom = domain.quantities['xmomentum'].explicit_update |
---|
54 | |
---|
55 | Elevation = domain.quantities['elevation'] |
---|
56 | Height = domain.quantities['height'] |
---|
57 | |
---|
58 | hc = Height.centroid_values |
---|
59 | b = Elevation.vertex_values |
---|
60 | |
---|
61 | x = domain.vertices |
---|
62 | g = domain.g |
---|
63 | |
---|
64 | x0 = x[:,0] |
---|
65 | x1 = x[:,1] |
---|
66 | |
---|
67 | b0 = b[:,0] |
---|
68 | b1 = b[:,1] |
---|
69 | |
---|
70 | bx = (b1-b0)/(x1-x0) |
---|
71 | |
---|
72 | xmom[:] = xmom - g*bx*hc |
---|
73 | |
---|
74 | def manning_friction(domain): |
---|
75 | """Apply (Manning) friction to water momentum |
---|
76 | """ |
---|
77 | |
---|
78 | from math import sqrt |
---|
79 | |
---|
80 | w = domain.quantities['stage'].centroid_values |
---|
81 | z = domain.quantities['elevation'].centroid_values |
---|
82 | h = w-z |
---|
83 | |
---|
84 | uh = domain.quantities['xmomentum'].centroid_values |
---|
85 | #vh = domain.quantities['ymomentum'].centroid_values |
---|
86 | eta = domain.quantities['friction'].centroid_values |
---|
87 | |
---|
88 | xmom_update = domain.quantities['xmomentum'].semi_implicit_update |
---|
89 | #ymom_update = domain.quantities['ymomentum'].semi_implicit_update |
---|
90 | |
---|
91 | N = domain.number_of_elements |
---|
92 | eps = domain.minimum_allowed_height |
---|
93 | g = domain.g |
---|
94 | |
---|
95 | for k in range(N): |
---|
96 | if eta[k] >= eps: |
---|
97 | if h[k] >= eps: |
---|
98 | #S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2)) |
---|
99 | S = -g * eta[k]**2 * uh[k] |
---|
100 | S /= h[k]**(7.0/3) |
---|
101 | |
---|
102 | #Update momentum |
---|
103 | xmom_update[k] += S*uh[k] |
---|
104 | #ymom_update[k] += S*vh[k] |
---|
105 | |
---|
106 | def linear_friction(domain): |
---|
107 | """Apply linear friction to water momentum |
---|
108 | |
---|
109 | Assumes quantity: 'linear_friction' to be present |
---|
110 | """ |
---|
111 | |
---|
112 | from math import sqrt |
---|
113 | |
---|
114 | w = domain.quantities['stage'].centroid_values |
---|
115 | z = domain.quantities['elevation'].centroid_values |
---|
116 | h = w-z |
---|
117 | |
---|
118 | uh = domain.quantities['xmomentum'].centroid_values |
---|
119 | # vh = domain.quantities['ymomentum'].centroid_values |
---|
120 | tau = domain.quantities['linear_friction'].centroid_values |
---|
121 | |
---|
122 | xmom_update = domain.quantities['xmomentum'].semi_implicit_update |
---|
123 | # ymom_update = domain.quantities['ymomentum'].semi_implicit_update |
---|
124 | |
---|
125 | N = domain.number_of_elements |
---|
126 | eps = domain.minimum_allowed_height |
---|
127 | g = domain.g #Not necessary? Why was this added? |
---|
128 | |
---|
129 | for k in range(N): |
---|
130 | if tau[k] >= eps: |
---|
131 | if h[k] >= eps: |
---|
132 | S = -tau[k]/h[k] |
---|
133 | |
---|
134 | #Update momentum |
---|
135 | xmom_update[k] += S*uh[k] |
---|
136 | # ymom_update[k] += S*vh[k] |
---|
137 | |
---|
138 | |
---|
139 | |
---|
140 | def check_forcefield(f): |
---|
141 | """Check that f is either |
---|
142 | 1: a callable object f(t,x,y), where x and y are vectors |
---|
143 | and that it returns an array or a list of same length |
---|
144 | as x and y |
---|
145 | 2: a scalar |
---|
146 | """ |
---|
147 | |
---|
148 | |
---|
149 | if callable(f): |
---|
150 | #N = 3 |
---|
151 | N = 2 |
---|
152 | #x = ones(3, numpy.float) |
---|
153 | #y = ones(3, numpy.float) |
---|
154 | x = ones(2, numpy.float) |
---|
155 | #y = ones(2, numpy.float) |
---|
156 | |
---|
157 | try: |
---|
158 | #q = f(1.0, x=x, y=y) |
---|
159 | q = f(1.0, x=x) |
---|
160 | except Exception, e: |
---|
161 | msg = 'Function %s could not be executed:\n%s' %(f, e) |
---|
162 | #FIXME: Reconsider this semantics |
---|
163 | raise msg |
---|
164 | |
---|
165 | try: |
---|
166 | q = numpy.array(q, numpy.float) |
---|
167 | except: |
---|
168 | msg = 'Return value from vector function %s could ' %f |
---|
169 | msg += 'not be converted into a numpy array of numpy.floats.\n' |
---|
170 | msg += 'Specified function should return either list or array.' |
---|
171 | raise msg |
---|
172 | |
---|
173 | #Is this really what we want? |
---|
174 | msg = 'Return vector from function %s ' %f |
---|
175 | msg += 'must have same lenght as input vectors' |
---|
176 | assert len(q) == N, msg |
---|
177 | |
---|
178 | else: |
---|
179 | try: |
---|
180 | f = float(f) |
---|
181 | except: |
---|
182 | msg = 'Force field %s must be either a scalar' %f |
---|
183 | msg += ' or a vector function' |
---|
184 | raise msg |
---|
185 | return f |
---|
186 | |
---|
187 | |
---|