 Timestamp:
 Feb 5, 2010, 5:24:39 PM (13 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_core/source/anuga/shallow_water_balanced/swb_domain.py
r7575 r7616 80 80 # set some defaults 81 81 # 82 self.set_timestepping_method( 1)82 self.set_timestepping_method(2) 83 83 self.set_default_order(2) 84 84 self.set_new_mannings_function(True) 85 85 self.set_centroid_transmissive_bc(True) 86 self.set_CFL(1.0) 87 self.set_beta(1.0) 86 88 87 89 ## … … 114 116 #Call correct module function (either from this module or Cextension) 115 117 118 #compute_fluxes(self) 119 120 #return 116 121 from swb_domain_ext import compute_fluxes_c 117 122 … … 129 134 self.flux_timestep = compute_fluxes_c(timestep, self, W, UH, VH, H, Z, U, V) 130 135 136 #print self.flux_timestep 137 131 138 ## 132 139 # @brief … … 150 157 """ 151 158 152 159 #Sww_domain.distribute_to_vertices_and_edges(self) 160 #return 161 153 162 #Shortcuts 154 163 W = self.quantities['stage'] … … 178 187 num.putmask(uh_C, h_C < 1.0e15, 0.0) 179 188 num.putmask(vh_C, h_C < 1.0e15, 0.0) 180 num.putmask(h_C, h_C < 1.0e15, 1.0e15) 181 182 u_C[:] = uh_C/h_C 183 v_C[:] = vh_C/h_C 184 185 for name in [ 'height', 'elevation', 'xvelocity', 'yvelocity' ]: 189 num.putmask(h_C, h_C < 1.0e15, 1.0e16) 190 191 H0 = 1.0e8 192 u_C[:] = uh_C/(h_C + H0/h_C) 193 v_C[:] = vh_C/(h_C + H0/h_C) 194 195 num.putmask(h_C, h_C < 1.0e15, 0.0) 196 197 for name in [ 'stage', 'height', 'xvelocity', 'yvelocity' ]: 186 198 Q = self.quantities[name] 187 199 if self._order_ == 1: … … 189 201 elif self._order_ == 2: 190 202 Q.extrapolate_second_order_and_limit_by_edge() 203 #Q.extrapolate_second_order_and_limit_by_vertex() 191 204 else: 192 205 raise 'Unknown order' … … 202 215 203 216 204 w_E[:] = z_E + h_E 205 206 #num.putmask(u_E, h_temp <= 0.0, 0.0) 207 #num.putmask(v_E, h_temp <= 0.0, 0.0) 208 #num.putmask(w_E, h_temp <= 0.0, z_E+h_E) 217 #minh_E = num.min(h_E) 218 #msg = 'min h_E = %g ' % minh_E 219 #assert minh_E >= 1.0e15, msg 220 221 z_E[:] = w_E  h_E 222 223 num.putmask(h_E, h_E <= 1.0e15, 0.0) 224 num.putmask(u_E, h_E <= 1.0e15, 0.0) 225 num.putmask(v_E, h_E <= 1.0e15, 0.0) 226 num.putmask(w_E, h_E <= 1.0e15, z_E) 209 227 #num.putmask(h_E, h_E <= 0.0, 0.0) 210 228 … … 246 264 247 265 266 267 268 ## 269 # @brief 270 def distribute_to_vertices_and_edges_h(self): 271 """Distribution from centroids to edges specific to the SWW eqn. 272 273 It will ensure that h (wz) is always nonnegative even in the 274 presence of steep bedslopes by taking a weighted average between shallow 275 and deep cases. 276 277 In addition, all conserved quantities get distributed as per either a 278 constant (order==1) or a piecewise linear function (order==2). 279 280 281 Precondition: 282 All conserved quantities defined at centroids and bed elevation defined at 283 edges. 284 285 Postcondition 286 Evolved quantities defined at vertices and edges 287 """ 288 289 290 #Shortcuts 291 W = self.quantities['stage'] 292 UH = self.quantities['xmomentum'] 293 VH = self.quantities['ymomentum'] 294 H = self.quantities['height'] 295 Z = self.quantities['elevation'] 296 U = self.quantities['xvelocity'] 297 V = self.quantities['yvelocity'] 298 299 #Arrays 300 w_C = W.centroid_values 301 uh_C = UH.centroid_values 302 vh_C = VH.centroid_values 303 z_C = Z.centroid_values 304 h_C = H.centroid_values 305 u_C = U.centroid_values 306 v_C = V.centroid_values 307 308 w_C[:] = num.maximum(w_C, z_C) 309 310 h_C[:] = w_C  z_C 311 312 313 assert num.min(h_C) >= 0 314 315 num.putmask(uh_C, h_C < 1.0e15, 0.0) 316 num.putmask(vh_C, h_C < 1.0e15, 0.0) 317 num.putmask(h_C, h_C < 1.0e15, 1.0e16) 318 319 u_C[:] = uh_C/h_C 320 v_C[:] = vh_C/h_C 321 322 num.putmask(h_C, h_C < 1.0e15, 0.0) 323 324 for name in [ 'stage', 'height', 'xvelocity', 'yvelocity' ]: 325 Q = self.quantities[name] 326 if self._order_ == 1: 327 Q.extrapolate_first_order() 328 elif self._order_ == 2: 329 Q.extrapolate_second_order_and_limit_by_edge() 330 #Q.extrapolate_second_order_and_limit_by_vertex() 331 else: 332 raise 'Unknown order' 333 334 335 w_E = W.edge_values 336 uh_E = UH.edge_values 337 vh_E = VH.edge_values 338 h_E = H.edge_values 339 z_E = Z.edge_values 340 u_E = U.edge_values 341 v_E = V.edge_values 342 343 344 #minh_E = num.min(h_E) 345 #msg = 'min h_E = %g ' % minh_E 346 #assert minh_E >= 1.0e15, msg 347 348 z_E[:] = w_E  h_E 349 350 num.putmask(h_E, h_E <= 1.0e8, 0.0) 351 num.putmask(u_E, h_E <= 1.0e8, 0.0) 352 num.putmask(v_E, h_E <= 1.0e8, 0.0) 353 num.putmask(w_E, h_E <= 1.0e8, z_E) 354 #num.putmask(h_E, h_E <= 0.0, 0.0) 355 356 uh_E[:] = u_E * h_E 357 vh_E[:] = v_E * h_E 358 359 """ 360 print '==========================================================' 361 print 'Time ', self.get_time() 362 print h_E 363 print uh_E 364 print vh_E 365 """ 366 367 # Compute vertex values by interpolation 368 for name in self.evolved_quantities: 369 Q = self.quantities[name] 370 Q.interpolate_from_edges_to_vertices() 371 372 373 w_V = W.vertex_values 374 uh_V = UH.vertex_values 375 vh_V = VH.vertex_values 376 z_V = Z.vertex_values 377 h_V = H.vertex_values 378 u_V = U.vertex_values 379 v_V = V.vertex_values 380 381 382 #w_V[:] = z_V + h_V 383 384 #num.putmask(u_V, h_V <= 0.0, 0.0) 385 #num.putmask(v_V, h_V <= 0.0, 0.0) 386 #num.putmask(w_V, h_V <= 0.0, z_V) 387 #num.putmask(h_V, h_V <= 0.0, 0.0) 388 389 uh_V[:] = u_V * h_V 390 vh_V[:] = v_V * h_V 391 392 393 394 248 395 ## 249 396 # @brief Code to let us use old shallow water domain BCs
Note: See TracChangeset
for help on using the changeset viewer.