source: trunk/anuga_work/development/2010-projects/anuga_1d/sww/sww_forcing_terms.py @ 7884

Last change on this file since 7884 was 7884, checked in by steve, 12 years ago

Moving 2010 project

File size: 4.7 KB
Line 
1#! /usr/bin/python
2
3
4"""
5Forcing 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
12def 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
49def 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
74def 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
106def 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
140def 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
Note: See TracBrowser for help on using the repository browser.