Changeset 221
- Timestamp:
- Aug 25, 2004, 5:26:37 PM (21 years ago)
- Location:
- inundation/ga/storm_surge/pyvolution
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/ga/storm_surge/pyvolution/config.py
r205 r221 59 59 #use_extensions = False #Do not use C-extensions 60 60 61 use_psyco = True #Use psyco optimisations62 #use_psyco = False #Do not use psyco optimisations61 #use_psyco = True #Use psyco optimisations 62 use_psyco = False #Do not use psyco optimisations 63 63 64 64 -
inundation/ga/storm_surge/pyvolution/quantity.py
r215 r221 78 78 v2 = self.vertex_values[i, 2] 79 79 80 self.centroid_values[i] = (v0 + v1 + v2)/3 80 self.centroid_values[i] = (v0 + v1 + v2)/3.0 81 81 82 82 self.edge_values[i, 0] = 0.5*(v1 + v2) … … 139 139 #Special case where we have only one neighbouring volume. 140 140 141 k0 = k #Self142 141 #Find index of the one neighbour 143 for k 1in self.domain.neighbours[k,:]:144 if k 1>= 0:142 for k0 in self.domain.neighbours[k,:]: 143 if k0 >= 0: 145 144 break 146 assert k1 != k0 147 assert k1 >= 0 145 146 assert k0 != k 147 assert k0 >= 0 148 149 k1 = k #Self 148 150 149 151 #Get data … … 180 182 x2, y2 = self.domain.centroids[k2] #V2 centroid 181 183 184 #if k < 4: 185 # print 'k, k0, k1, k2 = %d, %d, %d, %d' %(k, k0, k1, k2) 186 # #print 'k = %d: (%f, %f), (%f, %f), (%f, %f)'\ 187 # # %(k, x0, y0, x1, y1, x2, y2) 188 # print 'k = %d: (%f, %f, %f)' %(k, q0, q1, q2) 189 182 190 #Gradient 183 191 a[k], b[k] = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2) … … 266 274 267 275 a, b = self.compute_gradients() 276 ##for k in range(4): 277 ## print 'Gradients %d: %16.12f, %16.12f' %(k, a[k], b[k]) 268 278 269 279 V = self.domain.get_vertex_coordinates() -
inundation/ga/storm_surge/pyvolution/shallow_water.py
r215 r221 272 272 273 273 def distribute_to_vertices_and_edges(domain): 274 275 protect_against_infinitesimal_heights_centroid(domain)274 275 #print 'Calling distrib within SW' 276 276 if domain.order == 1: 277 protect_against_infinitesimal_heights_centroid(domain) 277 278 extrapolate_first_order(domain) 278 279 elif domain.order == 2: 280 protect_against_negative_heights_centroid(domain) 279 281 extrapolate_second_order(domain) 280 282 else: … … 291 293 minimum allowed height 292 294 """ 295 #FIXME: Used only in first order 293 296 294 297 #Water levels at centroids … … 316 319 317 320 321 def protect_against_negative_heights_centroid(domain): 322 """Adjust height and momentum at centroid if height is less than zero 323 """ 324 325 #FIXME: USed only in second order 326 #Water levels at centroids 327 wc = domain.quantities['level'].centroid_values 328 329 #Bed elevations at centroids 330 zc = domain.quantities['elevation'].centroid_values 331 332 #Water depths at centroids 333 hc = wc - zc 334 335 #Momentums at centroids 336 xmomc = domain.quantities['xmomentum'].centroid_values 337 ymomc = domain.quantities['ymomentum'].centroid_values 338 339 for k in range(domain.number_of_elements): 340 #Protect against infinitesimal heights and high velocities 341 if hc[k] < 0.0: 342 #Control level and height 343 wc[k] = zc[k]; hc[k] = 0.0 344 345 #Control momentum 346 xmomc[k] = ymomc[k] = 0.0 347 318 348 319 349 -
inundation/ga/storm_surge/pyvolution/test_shallow_water.py
r218 r221 960 960 #assert allclose(domain.max_timestep, 0.0396825396825) 961 961 #print domain.quantities['level'].centroid_values 962 963 964 def test_flatbed_first_order(self): 965 from mesh_factory import rectangular 966 from shallow_water import Domain,\ 967 Reflective_boundary, Dirichlet_boundary 968 969 from Numeric import array 970 971 #Create basic mesh 972 N = 8 973 points, vertices, boundary = rectangular(N, N) 974 975 #Create shallow water domain 976 domain = Domain(points, vertices, boundary) 977 domain.smooth = False 978 domain.visualise = False 979 domain.default_order=1 980 981 # Boundary conditions 982 Br = Reflective_boundary(domain) 983 Bd = Dirichlet_boundary([0.2,0.,0.]) 984 985 domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br}) 986 domain.check_integrity() 987 988 989 #Evolution 990 for t in domain.evolve(yieldstep = 0.02, finaltime = 0.5): 991 pass 992 #domain.write_time() 993 994 assert allclose(domain.min_timestep, 0.0140413643926) 995 assert allclose(domain.max_timestep, 0.0140947355753) 996 997 for i in range(3): 998 assert allclose(domain.quantities['level'].edge_values[:4,i], 999 [0.10730244,0.12337617,0.11200126,0.12605666]) 1000 1001 assert allclose(domain.quantities['xmomentum'].edge_values[:4,i], 1002 [0.07610894,0.06901572,0.07284461,0.06819712]) 1003 1004 assert allclose(domain.quantities['ymomentum'].edge_values[:4,i], 1005 [-0.0060238, -0.00157404, -0.00309633, -0.0001637]) 1006 1007 def test_flatbed_second_order(self): 1008 from mesh_factory import rectangular 1009 from shallow_water import Domain,\ 1010 Reflective_boundary, Dirichlet_boundary 1011 1012 from Numeric import array 1013 1014 #Create basic mesh 1015 N = 8 1016 points, vertices, boundary = rectangular(N, N) 1017 1018 #Create shallow water domain 1019 domain = Domain(points, vertices, boundary) 1020 domain.smooth = False 1021 domain.visualise = False 1022 domain.default_order=2 1023 1024 # Boundary conditions 1025 Br = Reflective_boundary(domain) 1026 Bd = Dirichlet_boundary([0.2,0.,0.]) 1027 1028 domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br}) 1029 domain.check_integrity() 1030 1031 #Evolution 1032 for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03): 1033 pass 1034 1035 1036 assert allclose(domain.min_timestep, 0.0210448446782) 1037 assert allclose(domain.max_timestep, 0.0210448446782) 1038 1039 #print domain.quantities['level'].vertex_values[:4,0] 1040 #print domain.quantities['xmomentum'].vertex_values[:4,0] 1041 #print domain.quantities['ymomentum'].vertex_values[:4,0] 1042 1043 assert allclose(domain.quantities['level'].vertex_values[:4,0], 1044 [0.00101913,0.05352143,0.00104852,0.05354394]) 1045 assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0], 1046 [0.00090581,0.03685719,0.00088303,0.03687313]) 1047 assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0], 1048 [-0.00139463,0.0006156,-0.00060364,0.00061827]) 1049 1050 1051 1052 def test_flatbed_second_order_distribute(self): 1053 #Use real data from pyvolution 2 1054 #painfully setup and extracted. 1055 from mesh_factory import rectangular 1056 from shallow_water import Domain,\ 1057 Reflective_boundary, Dirichlet_boundary 1058 1059 from Numeric import array 1060 1061 #Create basic mesh 1062 N = 8 1063 points, vertices, boundary = rectangular(N, N) 1064 1065 #Create shallow water domain 1066 domain = Domain(points, vertices, boundary) 1067 domain.smooth = False 1068 domain.visualise = False 1069 domain.default_order=domain.order=2 1070 1071 # Boundary conditions 1072 Br = Reflective_boundary(domain) 1073 Bd = Dirichlet_boundary([0.2,0.,0.]) 1074 1075 domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br}) 1076 domain.check_integrity() 1077 1078 1079 1080 for V in [False, True]: 1081 if V: 1082 #Set centroids as if system had been evolved 1083 L = zeros(2*N*N, Float) 1084 L[:32] = [7.21205592e-003, 5.35214298e-002, 1.00910824e-002, 1085 5.35439433e-002, 1.00910824e-002, 5.35439433e-002, 1086 1.00910824e-002, 5.35439433e-002, 1.00910824e-002, 1087 5.35439433e-002, 1.00910824e-002, 5.35439433e-002, 1088 1.00910824e-002, 5.35393928e-002, 1.02344264e-002, 1089 5.59605058e-002, 0.00000000e+000, 3.31027800e-004, 1090 0.00000000e+000, 4.37962142e-005, 0.00000000e+000, 1091 4.37962142e-005, 0.00000000e+000, 4.37962142e-005, 1092 0.00000000e+000, 4.37962142e-005, 0.00000000e+000, 1093 4.37962142e-005, 0.00000000e+000, 4.37962142e-005, 1094 0.00000000e+000, 5.57305948e-005] 1095 1096 X = zeros(2*N*N, Float) 1097 X[:32] = [6.48351607e-003, 3.68571894e-002, 8.50733285e-003, 1098 3.68731327e-002, 8.50733285e-003, 3.68731327e-002, 1099 8.50733285e-003, 3.68731327e-002, 8.50733285e-003, 1100 3.68731327e-002, 8.50733285e-003, 3.68731327e-002, 1101 8.50733285e-003, 3.68693861e-002, 8.65220973e-003, 1102 3.85055387e-002, 0.00000000e+000, 2.86060840e-004, 1103 0.00000000e+000, 3.58905503e-005, 0.00000000e+000, 1104 3.58905503e-005, 0.00000000e+000, 3.58905503e-005, 1105 0.00000000e+000, 3.58905503e-005, 0.00000000e+000, 1106 3.58905503e-005, 0.00000000e+000, 3.58905503e-005, 1107 0.00000000e+000, 4.57662812e-005] 1108 1109 Y = zeros(2*N*N, Float) 1110 Y[:32]=[-1.39463104e-003, 6.15600298e-004, -6.03637382e-004, 1111 6.18272251e-004, -6.03637382e-004, 6.18272251e-004, 1112 -6.03637382e-004, 6.18272251e-004, -6.03637382e-004, 1113 6.18272251e-004, -6.03637382e-004, 6.18272251e-004, 1114 -6.03637382e-004, 6.18599320e-004, -6.74622797e-004, 1115 -1.48934756e-004, 0.00000000e+000, -5.35079969e-005, 1116 0.00000000e+000, -2.57264987e-005, 0.00000000e+000, 1117 -2.57264987e-005, 0.00000000e+000, -2.57264987e-005, 1118 0.00000000e+000, -2.57264987e-005, 0.00000000e+000, 1119 -2.57264987e-005, 0.00000000e+000, -2.57264987e-005, 1120 0.00000000e+000, -2.57635178e-005] 1121 1122 1123 domain.set_quantity('level', L, 'centroids') 1124 domain.set_quantity('xmomentum', X, 'centroids') 1125 domain.set_quantity('ymomentum', Y, 'centroids') 1126 1127 domain.check_integrity() 1128 else: 1129 #Evolution 1130 for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03): 1131 pass 1132 assert allclose(domain.min_timestep, 0.0210448446782) 1133 assert allclose(domain.max_timestep, 0.0210448446782) 1134 1135 1136 #Centroids were correct but not vertices. 1137 #Hence the check of distribute below. 1138 assert allclose(domain.quantities['level'].centroid_values[:4], 1139 [0.00721206,0.05352143,0.01009108,0.05354394]) 1140 1141 assert allclose(domain.quantities['xmomentum'].centroid_values[:4], 1142 [0.00648352,0.03685719,0.00850733,0.03687313]) 1143 1144 assert allclose(domain.quantities['ymomentum'].centroid_values[:4], 1145 [-0.00139463,0.0006156,-0.00060364,0.00061827]) 1146 1147 #print 'C17=', domain.quantities['xmomentum'].centroid_values[17] 1148 #print 'C19=', domain.quantities['xmomentum'].centroid_values[19] 1149 1150 assert allclose(domain.quantities['xmomentum'].centroid_values[17], 1151 0.00028606084) 1152 1153 import copy 1154 XX = copy.copy(domain.quantities['xmomentum'].centroid_values) 1155 assert allclose(domain.quantities['xmomentum'].centroid_values, XX) 1156 1157 domain.distribute_to_vertices_and_edges() 1158 1159 assert allclose(domain.quantities['xmomentum'].centroid_values, XX) 1160 1161 assert allclose(domain.quantities['xmomentum'].centroid_values[17], 1162 0.00028606084) 1163 1164 1165 assert allclose(domain.quantities['level'].vertex_values[:4,0], 1166 [0.00101913,0.05352143,0.00104852,0.05354394]) 1167 1168 assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0], 1169 [-0.00139463,0.0006156,-0.00060364,0.00061827]) 1170 1171 #This was the culprit. First triangles vertex 0 had an 1172 #x-momentum of 0.0064835 instead of 0.00090581 and 1173 #third triangle had 0.00850733 instead of 0.00088303 1174 #print domain.quantities['xmomentum'].vertex_values[:4,0] 1175 1176 assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0], 1177 [0.00090581,0.03685719,0.00088303,0.03687313]) 1178 1179 962 1180 963 1181
Note: See TracChangeset
for help on using the changeset viewer.