Changeset 612
 Timestamp:
 Nov 22, 2004, 5:29:10 PM (19 years ago)
 Location:
 inundation/ga/storm_surge/pyvolution
 Files:

 3 edited
Legend:
 Unmodified
 Added
 Removed

inundation/ga/storm_surge/pyvolution/shallow_water.py
r610 r612 896 896 897 897 898 def check_forcefield(f): 899 """Check that f is either 900 1: a callable object of x,y,t, where x and y are vectors 901 and that it returns an array or a list of same length 902 as x and y 903 2: a scalar 904 """ 905 906 from Numeric import ones, Float, array 907 908 909 if callable(f): 910 N = 3 911 x = ones(3, Float) 912 y = ones(3, Float) 913 try: 914 q = f(x, y, 1.0) 915 except Exception, e: 916 msg = 'Function %s could not be executed:\n%s' %(f, e) 917 raise msg 918 919 try: 920 q = array(q).astype(Float) 921 except: 922 msg = 'Return value from vector function %s could ' %f 923 msg += 'not be converted into a Numeric array of floats.\n' 924 msg += 'Specified function should return either list or array.' 925 raise msg 926 927 msg = 'Return vector from function %s' %f 928 msg += 'must have same lenght as input vectors' 929 assert len(q) == N, msg 930 931 else: 932 try: 933 f = float(f) 934 except: 935 msg = 'Force field %s must be either a scalar' %f 936 msg += ' or a vector function' 937 raise msg 938 return f 939 898 940 899 941 class Wind_stress: … … 938 980 939 981 from config import rho_a, rho_w, eta_w 940 941 self.s = s 942 self.phi = phi 982 from Numeric import array, Float 983 984 self.speed = check_forcefield(s) 985 self.phi = check_forcefield(phi) 943 986 944 987 self.const = eta_w*rho_a/rho_w … … 950 993 951 994 from math import pi, cos, sin, sqrt 952 from Numeric import Float, ones 995 from Numeric import Float, ones, ArrayType 953 996 954 997 xmom_update = domain.quantities['xmomentum'].explicit_update … … 958 1001 t = domain.time 959 1002 960 if callable(self.s ):1003 if callable(self.speed): 961 1004 xc = domain.get_centroid_coordinates() 962 x = xc[:,0] 963 y = xc[:,1] 964 965 s_vec = self.s(x,y,t) 1005 s_vec = self.speed(xc[:,0], xc[:,1], t) 966 1006 else: 967 1007 #Assume s is a scalar 968 1008 969 1009 try: 970 s_vec = self.s * ones(N, Float)1010 s_vec = self.speed * ones(N, Float) 971 1011 except: 972 1012 msg = 'Speed must be either callable or a scalar: %s' %self.s … … 976 1016 if callable(self.phi): 977 1017 xc = domain.get_centroid_coordinates() 978 x = xc[:,0] 979 y = xc[:,1] 980 981 phi_vec = self.phi(x,y,t) 1018 phi_vec = self.phi(xc[:,0], xc[:,1], t) 982 1019 else: 983 1020 #Assume phi is a scalar … … 988 1025 msg = 'Angle must be either callable or a scalar: %s' %self.phi 989 1026 raise msg 1027 990 1028 991 1029 #This is the bit that should be written in C 
inundation/ga/storm_surge/pyvolution/test_shallow_water.py
r610 r612 30 30 31 31 return s 32 33 34 def scalar_func(x,y,t): 35 """Function that returns a scalar. 36 Used to test error message when Numeric array is expected 37 """ 38 39 return 17.7 32 40 33 41 … … 829 837 830 838 831 #FIXME: Do this when I get a minute of peace  goddammit 832 833 839 def test_wind_stress_error_condition(self): 840 """Test that windstress reacts properly when forcing functions 841 are wrong  e.g. returns a scalar 842 """ 843 844 from config import rho_a, rho_w, eta_w 845 from math import pi, cos, sin, sqrt 846 847 a = [0.0, 0.0] 848 b = [0.0, 2.0] 849 c = [2.0, 0.0] 850 d = [0.0, 4.0] 851 e = [2.0, 2.0] 852 f = [4.0, 0.0] 853 854 points = [a, b, c, d, e, f] 855 #bac, bce, ecf, dbe 856 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 857 858 domain = Domain(points, vertices) 859 860 #Flat surface with 1m of water 861 domain.set_quantity('elevation', 0) 862 domain.set_quantity('level', 1.0) 863 domain.set_quantity('friction', 0) 864 865 Br = Reflective_boundary(domain) 866 domain.set_boundary({'exterior': Br}) 867 868 869 domain.time = 5.54 #Take a random time (not zero) 870 871 #Setup only one forcing term, bad func 872 domain.forcing_terms = [] 873 874 try: 875 domain.forcing_terms.append(Wind_stress(s = scalar_func, 876 phi = angle)) 877 except AssertionError: 878 pass 879 else: 880 msg = 'Should have raised exception' 881 raise msg 882 883 884 try: 885 domain.forcing_terms.append(Wind_stress(s = speed, 886 phi = scalar_func)) 887 except AssertionError: 888 pass 889 else: 890 msg = 'Should have raised exception' 891 raise msg 892 893 try: 894 domain.forcing_terms.append(Wind_stress(s = speed, 895 phi = 'xx')) 896 except: 897 pass 898 else: 899 msg = 'Should have raised exception' 900 raise msg 834 901 835 902 
inundation/ga/storm_surge/pyvolution/wind_example_variable.py
r610 r612 13 13 #Create basic mesh 14 14 N = 80 15 len = 10016 points, vertices, boundary = rectangular(N, N, len , len)15 length = 100 16 points, vertices, boundary = rectangular(N, N, length, length) 17 17 18 18 #Create shallow water domain … … 35 35 Low speeds at center and edges 36 36 """ 37 38 from math import sqrt, exp, cos, pi 39 40 N = len(x) 41 s = 0*x #New array 42 43 c = (length/2, length/2) 37 44 38 from math import sqrt, exp, cos, pi 39 40 c = (len1/2, len2/2) 41 r = sqrt((x  c[0])**2 + (y  c[1])**2) 45 for k in range(N): 46 r = sqrt((x[k]  c[0])**2 + (y[k]  c[1])**2) 42 47 43 r /= max(len1,len2)48 r /= max(length,length) 44 49 45 factor = exp( (r0.15)**2 )50 factor = exp( (r0.15)**2 ) 46 51 47 #print x,y,factor52 s[k] = 4000 * factor * (cos(t*2*pi/150) + 2) 48 53 49 return 4000 * factor * (cos(t*2*pi/150) + 2)54 return s 50 55 51 56 … … 55 60 from math import sqrt, atan, pi 56 61 57 c = (len1/2, len2/2) 58 r = sqrt((x  c[0])**2 + (y  c[1])**2) 62 N = len(x) 63 a = 0*x #New array 64 65 c = (length/2, length/2) 66 for k in range(N): 67 r = sqrt((x[k]  c[0])**2 + (y[k]  c[1])**2) 68 69 xx = (x[k]  c[0])/length 70 yy = (y[k]  c[1])/length 71 72 angle = atan(yy/xx) 73 74 if xx < 0: 75 angle+=pi #atan in ]pi/2; pi/2[ 76 77 #Take normal direction 78 angle = pi/2 59 79 60 xx = (x  c[0])/len1 61 yy = (y  c[1])/len2 80 #Ensure positive radians 81 if angle < 0: 82 angle += 2*pi 62 83 63 angle = atan(yy/xx) 64 65 if xx < 0: 66 angle+=pi #atan in ]pi/2; pi/2[ 67 68 #Take normal direction 69 angle = pi/2 70 71 #Ensure positive radians 72 if angle < 0: 73 angle += 2*pi 74 75 return angle/pi*180 84 a[k] = angle/pi*180 85 86 return a 76 87 77 88 … … 82 93 def gust(x,y,t): 83 94 from math import sin, pi 95 from Numeric import zeros, ones, Float 96 97 N = len(x) 84 98 85 99 tt = sin(2*pi*t/200) 86 100 87 101 if tt > 0.9: 88 return 6000*tt 102 return 6000*tt*ones(N, Float) 89 103 else: 90 return 0104 return zeros(N, Float) 91 105 92 106 domain.forcing_terms.append(Wind_stress(gust, 25))
Note: See TracChangeset
for help on using the changeset viewer.