Changeset 5563 for anuga_work/development/anuga_1d/test_shallow_water.py
 Timestamp:
 Jul 23, 2008, 4:38:10 PM (15 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_work/development/anuga_1d/test_shallow_water.py
r5538 r5563 41 41 assert allclose( domain.quantities['stage'].explicit_update, stage_ud ) 42 42 assert allclose( domain.quantities['xmomentum'].explicit_update, xmom_ud ) 43 44 45 def test_local_flux_function(self): 46 normal = 1.0 47 ql = array([1.0, 2.0],Float) 48 qr = array([1.0, 2.0],Float) 49 zl = 0.0 50 zr = 0.0 51 52 edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr) 53 54 assert allclose(array([2.0, 8.9],Float), edgeflux) 55 assert allclose(5.1304951685, maxspeed) 56 57 normal = 1.0 58 ql = array([1.0, 2.0],Float) 59 qr = array([1.0, 2.0],Float) 60 zl = 0.0 61 zr = 0.0 62 63 edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr) 64 65 assert allclose(array([2.0, 8.9],Float), edgeflux) 66 assert allclose(5.1304951685, maxspeed) 67 68 def test_domain_flux_function(self): 69 normal = 1.0 70 ql = array([1.0, 2.0],Float) 71 qr = array([1.0, 2.0],Float) 72 zl = 0.0 73 zr = 0.0 74 75 edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr) 76 77 #print edgeflux 78 79 from shallow_water_domain import flux_function as domain_flux_function 80 81 domainedgeflux, domainmaxspeed = domain_flux_function(normal, ql,qr,zl,zr) 82 83 #print domainedgeflux 84 85 assert allclose(domainedgeflux, edgeflux) 86 assert allclose(domainmaxspeed, maxspeed) 87 43 88 44 89 … … 128 173 #m = domain.neighbour_edges[k,i] 129 174 m = domain.neighbour_vertices[k,i] 175 #print i, ' ' , m 130 176 #qr = [stage[n, m], xmom[n, m], ymom[n, m]] 131 177 qr[0] = stage[n, m] … … 169 215 170 216 217 def flux_function(normal, ql, qr, zl, zr): 218 """Compute fluxes between volumes for the shallow water wave equation 219 cast in terms of w = h+z using the 'central scheme' as described in 220 221 Kurganov, Noelle, Petrova. 'Semidiscrete CentralUpwind Schemes For 222 Hyperbolic Conservation Laws and HamiltonJacobi Equations'. 223 Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707740. 224 225 The implemented formula is given in equation (3.15) on page 714 226 227 Conserved quantities w, uh, are stored as elements 0 and 1 228 in the numerical vectors ql an qr. 229 230 Bed elevations zl and zr. 231 """ 232 233 from config import g, epsilon 234 from math import sqrt 235 from Numeric import array 236 237 #print 'ql',ql 238 239 #Align momentums with xaxis 240 #q_left = rotate(ql, normal, direction = 1) 241 #q_right = rotate(qr, normal, direction = 1) 242 q_left = ql 243 q_left[1] = q_left[1]*normal 244 q_right = qr 245 q_right[1] = q_right[1]*normal 246 247 #z = (zl+zr)/2 #Take average of field values 248 z = 0.5*(zl+zr) #Take average of field values 249 250 w_left = q_left[0] #w=h+z 251 h_left = w_leftz 252 uh_left = q_left[1] 253 254 if h_left < epsilon: 255 u_left = 0.0 #Could have been negative 256 h_left = 0.0 257 else: 258 u_left = uh_left/h_left 259 260 261 w_right = q_right[0] #w=h+z 262 h_right = w_rightz 263 uh_right = q_right[1] 264 265 266 if h_right < epsilon: 267 u_right = 0.0 #Could have been negative 268 h_right = 0.0 269 else: 270 u_right = uh_right/h_right 271 272 #vh_left = q_left[2] 273 #vh_right = q_right[2] 274 275 #print h_right 276 #print u_right 277 #print h_left 278 #print u_right 279 280 soundspeed_left = sqrt(g*h_left) 281 soundspeed_right = sqrt(g*h_right) 282 283 #Maximal wave speed 284 s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0) 285 286 #Minimal wave speed 287 s_min = min(u_leftsoundspeed_left, u_rightsoundspeed_right, 0) 288 289 #Flux computation 290 291 #flux_left = array([u_left*h_left, 292 # u_left*uh_left + 0.5*g*h_left**2]) 293 #flux_right = array([u_right*h_right, 294 # u_right*uh_right + 0.5*g*h_right**2]) 295 flux_left = array([u_left*h_left, 296 u_left*uh_left + 0.5*g*h_left*h_left]) 297 flux_right = array([u_right*h_right, 298 u_right*uh_right + 0.5*g*h_right*h_right]) 299 300 denom = s_maxs_min 301 if denom == 0.0: 302 edgeflux = array([0.0, 0.0]) 303 max_speed = 0.0 304 else: 305 edgeflux = (s_max*flux_left  s_min*flux_right)/denom 306 edgeflux += s_max*s_min*(q_rightq_left)/denom 307 308 edgeflux[1] = edgeflux[1]*normal 309 310 max_speed = max(abs(s_max), abs(s_min)) 311 312 return edgeflux, max_speed 313 171 314 172 315 #
Note: See TracChangeset
for help on using the changeset viewer.