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

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

Moving 2010 project

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.