source: anuga_work/development/sudi/sw_1d/periodic_waves/johns/util.py @ 7837

Last change on this file since 7837 was 7837, checked in by mungkasi, 14 years ago

Again, adding some codes for 1d problems on debris avalanche and periodic waves.

File size: 4.8 KB
Line 
1"""
2This module contains various auxiliary function used by pyvolution.
3"""
4
5def mean(x):
6    from Numeric import sum
7    return sum(x)/len(x)
8
9
10def gradient(x0, x1, q0, q1):
11
12    if q1-q0 != 0:
13        a = (q1-q0)/(x1-x0)
14    else:
15        a = 0
16       
17    return a
18
19def  minmod(beta_p,beta_m):
20    if (abs(beta_p) < abs(beta_m)) & (beta_p*beta_m > 0.0):
21        phi = beta_p
22    elif (abs(beta_m) < abs(beta_p)) & (beta_p*beta_m > 0.0):
23        phi = beta_m
24    else:
25        phi = 0.0
26    return phi
27
28def  minmod_kurganov(a,b,c):
29    from Numeric import sign
30    if sign(a)==sign(b)==sign(c):
31        return sign(a)*min(abs(a),abs(b),abs(c))
32    else:
33        return 0.0
34
35def  maxmod(a,b):
36    if (abs(a) > abs(b)) & (a*b > 0.0):
37        phi = a
38    elif (abs(b) > abs(a)) & (a*b > 0.0):
39        phi = b
40    else:
41        phi = 0.0
42    return phi
43
44def vanleer(a,b):
45    if abs(a)+abs(b) > 1e-12:
46        return (a*abs(b)+abs(a)*b)/(abs(a)+abs(b))
47    else:
48        return 0.0
49
50def vanalbada(a,b):
51    if a*a+b*b > 1e-12:
52        return (a*a*b+a*b*b)/(a*a+b*b)
53    else:
54        return 0.0
55
56def calculate_wetted_area(x1,x2,z1,z2,w1,w2):
57    if (w1 > z1) & (w2 < z2) & (z1 <= z2):
58        x = ((w2-z1)*(x2-x1)+x1*(z2-z1)-x2*(w2-w1))/(z2-z1+w1-w2)
59        A = 0.5*(w1-z1)*(x-x1)
60        L = x-x1
61    elif (w1 < z1) & (w2 > z2) & (z1 < z2):
62        x = ((w2-z1)*(x2-x1)+x1*(z2-z1)-x2*(w2-w1))/(z2-z1+w1-w2)
63        A = 0.5*(w2-z2)*(x2-x)
64        L = x2-x
65    elif (w1 < z1) & (w2 > z2) & (z1 >= z2):
66        x = ((w1-z2)*(x2-x1)+x2*(z2-z1)-x1*(w2-w1))/(z2-z1+w1-w2)
67        A = 0.5*(w2-z2)*(x2-x)
68        L = x2-x
69    elif (w1 > z1) & (w2 < z2) & (z1 > z2):
70        x = ((w1-z2)*(x2-x1)+x2*(z2-z1)-x1*(w2-w1))/(z2-z1+w1-w2)
71        A = 0.5*(w1-z1)*(x-x1)
72        L = x-x1
73    elif (w1 <= z1) & (w2 <= z2):
74        A = 0.0
75    elif (w1 == z1) & (w2 > z2) & (z2 < z1):
76        A = 0.5*(x2-x1)*(w2-z2)
77    elif (w2 == z2) & (w1 > z1) & (z1 < z2):
78        A = 0.5*(x2-x1)*(w1-z1)
79    return A
80
81
82def calculate_new_wet_area(x1,x2,z1,z2,A):
83    from Numeric import sqrt
84    min_centroid_height = 1.0e-3
85    # Assumes reconstructed stage flat in a wetted cell
86    M = (z2-z1)/(x2-x1)
87    L = (x2-x1)
88    min_area = min_centroid_height*L
89    max_area = 0.5*(x2-x1)*abs(z2-z1)
90    if A < max_area:
91        if (z1 < z2):
92            x = sqrt(2*A/M)+x1
93            wet_len = x-x1
94            wc = z1 + sqrt(M*2*A)
95        elif (z2 < z1):
96            x = -sqrt(-2*A/M)+x2
97            wet_len = x2-x
98            wc = z2+sqrt(-M*2*A)
99        else:
100            wc = A/L+0.5*(z1+z2)
101            wet_len = x2-x1
102    else:
103        wc = 0.5*(z1+z2)+A/L
104        wet_len = x2-x1
105           
106    return wc,wet_len
107
108def calculate_new_wet_area_analytic(x1,x2,z1,z2,A,t):
109    min_centroid_height = 1.0e-3
110    # Assumes reconstructed stage flat in a wetted cell
111    M = (z2-z1)/(x2-x1)
112    L = (x2-x1)
113    min_area = min_centroid_height*L
114    max_area = 0.5*(x2-x1)*abs(z2-z1)
115    w1,uh1 = analytic_cannal(x1,t)
116    w2,uh2 = analytic_cannal(x2,t)
117    if (w1 > z1) & (w2 < z2) & (z1 <= z2):
118        print "test1"
119        x = ((w2-z1)*(x2-x1)+x1*(z2-z1)-x2*(w2-w1))/(z2-z1+w1-w2)
120        wet_len = x-x1
121    elif (w1 < z1) & (w2 > z2) & (z1 < z2):
122        print "test2"
123        x = ((w2-z1)*(x2-x1)+x1*(z2-z1)-x2*(w2-w1))/(z2-z1+w1-w2)
124        wet_len = x2-x
125    elif (w1 < z1) & (w2 > z2) & (z1 >= z2):
126        print "test3"
127        x = ((w1-z2)*(x2-x1)+x2*(z2-z1)-x1*(w2-w1))/(z2-z1+w1-w2)
128        wet_len = x2-x
129    elif (w1 > z1) & (w2 < z2) & (z1 > z2):
130        print "test4"
131        x = ((w1-z2)*(x2-x1)+x2*(z2-z1)-x1*(w2-w1))/(z2-z1+w1-w2)
132        wet_len = x-x1
133    elif (w1 >= z1) & (w2 >= z2):
134        print "test5"
135        wet_len = x2-x1
136    else: #(w1 <= z1) & (w2 <= z2)
137        print "test5"
138        if (w1 > z1) | (w2 > z2):
139            print "ERROR"
140        wet_len = x2-x1       
141    return w1,w2,wet_len,uh1,uh2
142
143def analytic_cannal(C,t):
144    from Numeric import zeros, Float,sqrt,sin,cos
145
146   
147    #N = len(C)
148    #u = zeros(N,Float)    ## water velocity
149    #h = zeros(N,Float)    ## water depth
150    x = C
151    g = 9.81
152
153
154    ## Define Basin Bathymetry
155    #z_b = zeros(N,Float) ## elevation of basin
156    #z = zeros(N,Float)   ## elevation of water surface
157    z_infty = 10.0       ## max equilibrium water depth at lowest point.
158    L_x = 2500.0         ## width of channel
159
160    A0 = 0.5*L_x                  ## determines amplitudes of oscillations
161    omega = sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation
162
163    x1 = A0*cos(omega*t)-L_x # left shoreline
164    x2 = A0*cos(omega*t)+L_x # right shoreline
165    if (x >=x1) & (x <= x2):
166        z_b = z_infty*(x**2/L_x**2) ## or A0*cos(omega*t)\pmL_x
167        u = -A0*omega*sin(omega*t)
168        z = z_infty+2*A0*z_infty/L_x*cos(omega*t)*(x/L_x-0.5*A0/(L_x)*cos(omega*t))
169    else:
170       z_b = z_infty*(x**2/L_x**2)
171       u=0.0
172       z = z_b
173    h = z-z_b
174    return z,u*h
Note: See TracBrowser for help on using the repository browser.