source: anuga_work/development/2010-projects/anuga_1d/sww/wet_dry.py @ 7839

Last change on this file since 7839 was 7839, checked in by steve, 14 years ago

Changing name of 1d projects so that it will be easy to moveto the trunk

File size: 3.8 KB
Line 
1"""
2This module contains various auxiliary function used
3to ascertain wet dry areas.
4"""
5
6
7import numpy
8
9def calculate_wetted_area(x1,x2,z1,z2,w1,w2):
10    if (w1 > z1) & (w2 < z2) & (z1 <= z2):
11        x = ((w2-z1)*(x2-x1)+x1*(z2-z1)-x2*(w2-w1))/(z2-z1+w1-w2)
12        A = 0.5*(w1-z1)*(x-x1)
13        L = x-x1
14    elif (w1 < z1) & (w2 > z2) & (z1 < z2):
15        x = ((w2-z1)*(x2-x1)+x1*(z2-z1)-x2*(w2-w1))/(z2-z1+w1-w2)
16        A = 0.5*(w2-z2)*(x2-x)
17        L = x2-x
18    elif (w1 < z1) & (w2 > z2) & (z1 >= z2):
19        x = ((w1-z2)*(x2-x1)+x2*(z2-z1)-x1*(w2-w1))/(z2-z1+w1-w2)
20        A = 0.5*(w2-z2)*(x2-x)
21        L = x2-x
22    elif (w1 > z1) & (w2 < z2) & (z1 > z2):
23        x = ((w1-z2)*(x2-x1)+x2*(z2-z1)-x1*(w2-w1))/(z2-z1+w1-w2)
24        A = 0.5*(w1-z1)*(x-x1)
25        L = x-x1
26    elif (w1 <= z1) & (w2 <= z2):
27        A = 0.0
28    elif (w1 == z1) & (w2 > z2) & (z2 < z1):
29        A = 0.5*(x2-x1)*(w2-z2)
30    elif (w2 == z2) & (w1 > z1) & (z1 < z2):
31        A = 0.5*(x2-x1)*(w1-z1)
32    return A
33
34
35def calculate_new_wet_area(x1,x2,z1,z2,A):
36    from numpy import sqrt
37    min_centroid_height = 1.0e-3
38    # Assumes reconstructed stage flat in a wetted cell
39    M = (z2-z1)/(x2-x1)
40    L = (x2-x1)
41    min_area = min_centroid_height*L
42    max_area = 0.5*(x2-x1)*abs(z2-z1)
43    if A < max_area:
44        if (z1 < z2):
45            x = sqrt(2*A/M)+x1
46            wet_len = x-x1
47            wc = z1 + sqrt(M*2*A)
48        elif (z2 < z1):
49            x = -sqrt(-2*A/M)+x2
50            wet_len = x2-x
51            wc = z2+sqrt(-M*2*A)
52        else:
53            wc = A/L+0.5*(z1+z2)
54            wet_len = x2-x1
55    else:
56        wc = 0.5*(z1+z2)+A/L
57        wet_len = x2-x1
58           
59    return wc,wet_len
60
61def calculate_new_wet_area_analytic(x1,x2,z1,z2,A,t):
62    min_centroid_height = 1.0e-3
63    # Assumes reconstructed stage flat in a wetted cell
64    M = (z2-z1)/(x2-x1)
65    L = (x2-x1)
66    min_area = min_centroid_height*L
67    max_area = 0.5*(x2-x1)*abs(z2-z1)
68    w1,uh1 = analytic_cannal(x1,t)
69    w2,uh2 = analytic_cannal(x2,t)
70    if (w1 > z1) & (w2 < z2) & (z1 <= z2):
71        print "test1"
72        x = ((w2-z1)*(x2-x1)+x1*(z2-z1)-x2*(w2-w1))/(z2-z1+w1-w2)
73        wet_len = x-x1
74    elif (w1 < z1) & (w2 > z2) & (z1 < z2):
75        print "test2"
76        x = ((w2-z1)*(x2-x1)+x1*(z2-z1)-x2*(w2-w1))/(z2-z1+w1-w2)
77        wet_len = x2-x
78    elif (w1 < z1) & (w2 > z2) & (z1 >= z2):
79        print "test3"
80        x = ((w1-z2)*(x2-x1)+x2*(z2-z1)-x1*(w2-w1))/(z2-z1+w1-w2)
81        wet_len = x2-x
82    elif (w1 > z1) & (w2 < z2) & (z1 > z2):
83        print "test4"
84        x = ((w1-z2)*(x2-x1)+x2*(z2-z1)-x1*(w2-w1))/(z2-z1+w1-w2)
85        wet_len = x-x1
86    elif (w1 >= z1) & (w2 >= z2):
87        print "test5"
88        wet_len = x2-x1
89    else: #(w1 <= z1) & (w2 <= z2)
90        print "test5"
91        if (w1 > z1) | (w2 > z2):
92            print "ERROR"
93        wet_len = x2-x1       
94    return w1,w2,wet_len,uh1,uh2
95
96def analytic_cannal(C,t):
97
98    import math
99    #N = len(C)
100    #u = zeros(N,numpy.float)    ## water velocity
101    #h = zeros(N,numpy.float)    ## water depth
102    x = C
103    g = 9.81
104
105
106    ## Define Basin Bathymetry
107    #z_b = zeros(N,numpy.float) ## elevation of basin
108    #z = zeros(N,numpy.float)   ## elevation of water surface
109    z_infty = 10.0       ## max equilibrium water depth at lowest point.
110    L_x = 2500.0         ## width of channel
111
112    A0 = 0.5*L_x                  ## determines amplitudes of oscillations
113    omega = math.sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation
114
115    x1 = A0*cos(omega*t)-L_x # left shoreline
116    x2 = A0*cos(omega*t)+L_x # right shoreline
117    if (x >=x1) & (x <= x2):
118        z_b = z_infty*(x**2/L_x**2) ## or A0*cos(omega*t)\pmL_x
119        u = -A0*omega*sin(omega*t)
120        z = z_infty+2*A0*z_infty/L_x*cos(omega*t)*(x/L_x-0.5*A0/(L_x)*cos(omega*t))
121    else:
122       z_b = z_infty*(x**2/L_x**2)
123       u=0.0
124       z = z_b
125    h = z-z_b
126    return z,u*h
Note: See TracBrowser for help on using the repository browser.