Changeset 7835 for anuga_work/development/flow_1d/generic_1d/util.py
 Timestamp:
 Jun 14, 2010, 9:50:26 PM (13 years ago)
 File:

 1 moved
Legend:
 Unmodified
 Added
 Removed

anuga_work/development/flow_1d/generic_1d/util.py
r7834 r7835 56 56 return 0.0 57 57 58 def calculate_wetted_area(x1,x2,z1,z2,w1,w2):59 if (w1 > z1) & (w2 < z2) & (z1 <= z2):60 x = ((w2z1)*(x2x1)+x1*(z2z1)x2*(w2w1))/(z2z1+w1w2)61 A = 0.5*(w1z1)*(xx1)62 L = xx163 elif (w1 < z1) & (w2 > z2) & (z1 < z2):64 x = ((w2z1)*(x2x1)+x1*(z2z1)x2*(w2w1))/(z2z1+w1w2)65 A = 0.5*(w2z2)*(x2x)66 L = x2x67 elif (w1 < z1) & (w2 > z2) & (z1 >= z2):68 x = ((w1z2)*(x2x1)+x2*(z2z1)x1*(w2w1))/(z2z1+w1w2)69 A = 0.5*(w2z2)*(x2x)70 L = x2x71 elif (w1 > z1) & (w2 < z2) & (z1 > z2):72 x = ((w1z2)*(x2x1)+x2*(z2z1)x1*(w2w1))/(z2z1+w1w2)73 A = 0.5*(w1z1)*(xx1)74 L = xx175 elif (w1 <= z1) & (w2 <= z2):76 A = 0.077 elif (w1 == z1) & (w2 > z2) & (z2 < z1):78 A = 0.5*(x2x1)*(w2z2)79 elif (w2 == z2) & (w1 > z1) & (z1 < z2):80 A = 0.5*(x2x1)*(w1z1)81 return A82 83 84 def calculate_new_wet_area(x1,x2,z1,z2,A):85 from numpy import sqrt86 min_centroid_height = 1.0e387 # Assumes reconstructed stage flat in a wetted cell88 M = (z2z1)/(x2x1)89 L = (x2x1)90 min_area = min_centroid_height*L91 max_area = 0.5*(x2x1)*abs(z2z1)92 if A < max_area:93 if (z1 < z2):94 x = sqrt(2*A/M)+x195 wet_len = xx196 wc = z1 + sqrt(M*2*A)97 elif (z2 < z1):98 x = sqrt(2*A/M)+x299 wet_len = x2x100 wc = z2+sqrt(M*2*A)101 else:102 wc = A/L+0.5*(z1+z2)103 wet_len = x2x1104 else:105 wc = 0.5*(z1+z2)+A/L106 wet_len = x2x1107 108 return wc,wet_len109 110 def calculate_new_wet_area_analytic(x1,x2,z1,z2,A,t):111 min_centroid_height = 1.0e3112 # Assumes reconstructed stage flat in a wetted cell113 M = (z2z1)/(x2x1)114 L = (x2x1)115 min_area = min_centroid_height*L116 max_area = 0.5*(x2x1)*abs(z2z1)117 w1,uh1 = analytic_cannal(x1,t)118 w2,uh2 = analytic_cannal(x2,t)119 if (w1 > z1) & (w2 < z2) & (z1 <= z2):120 print "test1"121 x = ((w2z1)*(x2x1)+x1*(z2z1)x2*(w2w1))/(z2z1+w1w2)122 wet_len = xx1123 elif (w1 < z1) & (w2 > z2) & (z1 < z2):124 print "test2"125 x = ((w2z1)*(x2x1)+x1*(z2z1)x2*(w2w1))/(z2z1+w1w2)126 wet_len = x2x127 elif (w1 < z1) & (w2 > z2) & (z1 >= z2):128 print "test3"129 x = ((w1z2)*(x2x1)+x2*(z2z1)x1*(w2w1))/(z2z1+w1w2)130 wet_len = x2x131 elif (w1 > z1) & (w2 < z2) & (z1 > z2):132 print "test4"133 x = ((w1z2)*(x2x1)+x2*(z2z1)x1*(w2w1))/(z2z1+w1w2)134 wet_len = xx1135 elif (w1 >= z1) & (w2 >= z2):136 print "test5"137 wet_len = x2x1138 else: #(w1 <= z1) & (w2 <= z2)139 print "test5"140 if (w1 > z1)  (w2 > z2):141 print "ERROR"142 wet_len = x2x1143 return w1,w2,wet_len,uh1,uh2144 145 def analytic_cannal(C,t):146 147 import math148 #N = len(C)149 #u = zeros(N,numpy.float) ## water velocity150 #h = zeros(N,numpy.float) ## water depth151 x = C152 g = 9.81153 154 155 ## Define Basin Bathymetry156 #z_b = zeros(N,numpy.float) ## elevation of basin157 #z = zeros(N,numpy.float) ## elevation of water surface158 z_infty = 10.0 ## max equilibrium water depth at lowest point.159 L_x = 2500.0 ## width of channel160 161 A0 = 0.5*L_x ## determines amplitudes of oscillations162 omega = math.sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation163 164 x1 = A0*cos(omega*t)L_x # left shoreline165 x2 = A0*cos(omega*t)+L_x # right shoreline166 if (x >=x1) & (x <= x2):167 z_b = z_infty*(x**2/L_x**2) ## or A0*cos(omega*t)\pmL_x168 u = A0*omega*sin(omega*t)169 z = z_infty+2*A0*z_infty/L_x*cos(omega*t)*(x/L_x0.5*A0/(L_x)*cos(omega*t))170 else:171 z_b = z_infty*(x**2/L_x**2)172 u=0.0173 z = z_b174 h = zz_b175 return z,u*h
Note: See TracChangeset
for help on using the changeset viewer.