- Timestamp:
- Jun 15, 2010, 5:34:24 PM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/development/2010-projects/anuga_1d/channel/channel_domain.py
r7839 r7842 3 3 the shallow water wave equation. 4 4 5 This module contains a specialisation of class Domain from module domain.py 6 consisting of methods specific to the Shallow Water Wave Equation 5 This module contains a specialisation of class Generic_domain from module 6 generic_domain.py 7 consisting of methods specific to Channel flow using the Shallow Water Wave Equation 7 8 8 9 This particular modification of the Domain class implements the ability to … … 14 15 15 16 where 16 ------------!!!! NOTE THIS NEEDS UPDATING !!!!------------------ 17 U = [A, Q] 18 E = [Q, Q^2/A + g h^2/2]17 18 U = [A, Q] = [b*h, u*b*h] 19 E = [Q, Q^2/A + g*b*h^2/2] 19 20 S represents source terms forcing the system 20 (e.g. gravity, friction, wind stress, ...) 21 (e.g. gravity, boundary_stree, friction, wind stress, ...) 22 gravity = -g*b*h*z_x 23 boundary_stress = 1/2*g*b_x*h^2 21 24 22 25 and _t, _x, _y denote the derivative with respect to t, x and y respectiely. … … 25 28 26 29 symbol variable name explanation 30 A area Wetted area = b*h 31 Q discharge flux of water = u*b*h 27 32 x x horizontal distance from origin [m] 28 33 z elevation elevation of bed on which flow is modelled [m] … … 31 36 u speed in the x direction [m/s] 32 37 uh xmomentum momentum in the x direction [m^2/s] 33 38 b width width of channel 34 39 eta mannings friction coefficient [to appear] 35 40 nu wind stress coefficient [to appear] 36 41 37 The conserved quantities are w, uh42 The conserved quantities are A, Q 38 43 -------------------------------------------------------------------------- 39 44 For details see e.g. … … 48 53 49 54 50 from domain import * 51 Generic_Domain = Domain #Rename 55 from anuga_1d.generic.generic_domain import * 56 import numpy 57 52 58 53 59 #Shallow water domain 54 class Domain(Generic_ Domain):60 class Domain(Generic_domain): 55 61 56 62 def __init__(self, coordinates, boundary = None, tagged_elements = None): … … 59 65 evolved_quantities = ['area', 'discharge', 'elevation', 'height', 'velocity','width','stage'] 60 66 other_quantities = ['friction'] 61 Generic_ Domain.__init__(self,67 Generic_domain.__init__(self, 62 68 coordinates = coordinates, 63 69 boundary = boundary, … … 67 73 tagged_elements = tagged_elements) 68 74 69 from config import minimum_allowed_height, g, h075 from anuga_1d.config import minimum_allowed_height, g, h0 70 76 self.minimum_allowed_height = minimum_allowed_height 71 77 self.g = g … … 74 80 self.discontinousb = False 75 81 76 82 83 # forcing terms gravity and boundary stress are included in the flux calculation 77 84 #self.forcing_terms.append(gravity) 78 85 #self.forcing_terms.append(boundary_stress) … … 88 95 89 96 #Reduction operation for get_vertex_values 90 from util import mean97 from anuga_1d.generic.util import mean 91 98 self.reduction = mean 92 99 #self.reduction = min #Looks better near steep slopes … … 94 101 self.set_quantities_to_be_stored(['area','discharge']) 95 102 96 self.__doc__ = 'channel_domain _Ab'103 self.__doc__ = 'channel_domain' 97 104 98 105 self.check_integrity() … … 118 125 msg = 'Fifth evolved quantity must be "velocity"' 119 126 assert self.evolved_quantities[4] == 'velocity', msg 120 msg = ' Fifth evolved quantity must be "width"'127 msg = 'Sixth evolved quantity must be "width"' 121 128 assert self.evolved_quantities[5] == 'width', msg 122 123 Generic_Domain.check_integrity(self) 129 msg = 'Seventh evolved quantity must be "stage"' 130 assert self.evolved_quantities[6] == 'stage', msg 131 132 Generic_domain.check_integrity(self) 124 133 125 134 def compute_fluxes(self): … … 140 149 #----------------------------------- 141 150 def compute_fluxes_channel(domain): 142 from Numeric import zeros, Float143 151 import sys 144 145 146 152 timestep = float(sys.maxint) 147 153 148 area = domain.quantities['area']149 discharge 150 bed = domain.quantities['elevation']151 height = domain.quantities['height']152 velocity = domain.quantities['velocity']153 width = domain.quantities['width']154 155 156 from channel_domain_ext import compute_fluxes_channel_ext154 area = domain.quantities['area'] 155 discharge = domain.quantities['discharge'] 156 bed = domain.quantities['elevation'] 157 height = domain.quantities['height'] 158 velocity = domain.quantities['velocity'] 159 width = domain.quantities['width'] 160 161 162 from anuga_1d.channel.channel_domain_ext import compute_fluxes_channel_ext 157 163 domain.flux_timestep = compute_fluxes_channel_ext(timestep,domain,area,discharge,bed,height,velocity,width) 158 164 … … 168 174 169 175 import sys 170 from Numeric import zeros, Float 171 from config import epsilon, h0 172 ## linearb(domain) 173 174 176 from anuga_1d.config import epsilon, h0 175 177 176 178 … … 178 180 179 181 #Shortcuts 180 Area= domain.quantities['area']181 Discharge = domain.quantities['discharge']182 Bed= domain.quantities['elevation']183 Height= domain.quantities['height']184 Velocity= domain.quantities['velocity']185 Width= domain.quantities['width']186 Stage= domain.quantities['stage']182 area = domain.quantities['area'] 183 discharge = domain.quantities['discharge'] 184 bed = domain.quantities['elevation'] 185 height = domain.quantities['height'] 186 velocity = domain.quantities['velocity'] 187 width = domain.quantities['width'] 188 stage = domain.quantities['stage'] 187 189 188 190 #Arrays 189 a_C = Area.centroid_values190 d_C = Discharge.centroid_values191 z_C = Bed.centroid_values192 h_C = Height.centroid_values193 u_C = Velocity.centroid_values194 b_C = Width.centroid_values195 w_C = Stage.centroid_values191 a_C = area.centroid_values 192 d_C = discharge.centroid_values 193 z_C = bed.centroid_values 194 h_C = height.centroid_values 195 u_C = velocity.centroid_values 196 b_C = width.centroid_values 197 w_C = stage.centroid_values 196 198 197 199 if domain.setstageflag: 198 for i in range(len(a_C)): 199 a_C[i]=(w_C[i]-z_C[i])*b_C[i] 200 200 a_C[:,] = (w_C-z_C)*b_C 201 201 domain.setstageflag = False 202 202 … … 205 205 206 206 207 h0 = 1.0e-12 208 #print id(h_C) 209 for i in range(N): 210 211 if a_C[i] <= h0: 212 a_C[i] = 0.0 213 h_C[i] = 0.0 214 d_C[i] = 0.0 215 u_C[i] = 0.0 216 w_C[i] = z_C[i] 217 218 219 220 else: 221 222 if b_C[i]<=h0: 223 a_C[i] = 0.0 224 h_C[i] = 0.0 225 d_C[i] = 0.0 226 u_C[i] = 0.0 227 w_C[i] = z_C[i] 228 229 else: 230 h_C[i] = a_C[i]/(b_C[i]+h0/b_C[i]) 231 w_C[i] = h_C[i]+z_C[i] 232 u_C[i] = d_C[i]/(a_C[i]+h0/a_C[i]) 207 # FIXME (SGR): Replace with numpy.where to speed up code 208 h0 = 1.0e-12 209 210 h_C[:] = numpy.where( (a_C>h0) | (b_C>h0), a_C/(b_C+h0/b_C), 0.0 ) 211 u_C[:] = numpy.where( (a_C>h0) | (b_C>h0), d_C/(a_C+h0/a_C), 0.0 ) 212 213 a_C[:] = numpy.where( (a_C>h0) | (b_C>h0), a_C, 0.0 ) 214 b_C[:] = numpy.where( (a_C>h0) | (b_C>h0), b_C, 0.0 ) 215 d_C[:] = numpy.where( (a_C>h0) | (b_C>h0), d_C, 0.0 ) 216 217 w_C[:] = h_C + z_C 218 219 # for i in range(N): 220 # 221 # if a_C[i] <= h0: 222 # a_C[i] = 0.0 223 # h_C[i] = 0.0 224 # d_C[i] = 0.0 225 # u_C[i] = 0.0 226 # w_C[i] = z_C[i] 227 # else: 228 # 229 # if b_C[i]<=h0: 230 # a_C[i] = 0.0 231 # h_C[i] = 0.0 232 # d_C[i] = 0.0 233 # u_C[i] = 0.0 234 # w_C[i] = z_C[i] 235 # 236 # else: 237 # h_C[i] = a_C[i]/(b_C[i]+h0/b_C[i]) 238 # w_C[i] = h_C[i]+z_C[i] 239 # u_C[i] = d_C[i]/(a_C[i]+h0/a_C[i]) 233 240 234 235 236 241 242 bed.extrapolate_second_order() 243 237 244 for name in ['velocity','stage']: 238 245 Q = domain.quantities[name] … … 243 250 else: 244 251 raise 'Unknown order' 245 a_V = domain.quantities['area'].vertex_values 246 w_V = domain.quantities['stage'].vertex_values 247 z_V = domain.quantities['elevation'].vertex_values 248 h_V = domain.quantities['height'].vertex_values 249 u_V = domain.quantities['velocity'].vertex_values 250 d_V = domain.quantities['discharge'].vertex_values 251 b_V = domain.quantities['width'].vertex_values 252 253 254 a_V = area.vertex_values 255 w_V = stage.vertex_values 256 z_V = bed.vertex_values 257 h_V = height.vertex_values 258 u_V = velocity.vertex_values 259 d_V = discharge.vertex_values 260 b_V = width.vertex_values 252 261 253 262 254 255 for i in range(len(h_C)): 256 for j in range(2): 257 ## if b_V[i,j] < h0 : 258 ## a_V[i,j]=0 259 ## h_V[i,j]=0 260 ## d_V[i,j]=0 261 ## u_V[i,j]=0 262 ## else: 263 264 h_V[i,j] = w_V[i,j]-z_V[i,j] 265 if h_V[i,j]<h0: 266 h_V[i,j]=0 267 w_V[i,j]=z_V[i,j] 268 a_V[i,j] = b_V[i,j]*h_V[i,j] 269 d_V[i,j]=u_V[i,j]*a_V[i,j] 270 271 272 273 274 263 #FIXME (SGR): replace with numpy.where 264 265 266 h_V[:] = w_V-z_V 267 268 w_V[:] = numpy.where(h_V > h0, w_V, z_V) 269 h_V[:] = numpy.where(h_V > h0, h_V, 0.0) 270 a_V[:] = b_V*h_V 271 d_V[:] = u_V*a_V 272 273 # for i in range(len(h_C)): 274 # for j in range(2): 275 # h_V[i,j] = w_V[i,j]-z_V[i,j] 276 # if h_V[i,j]<h0: 277 # h_V[i,j]=0 278 # w_V[i,j]=z_V[i,j] 279 # a_V[i,j] = b_V[i,j]*h_V[i,j] 280 # d_V[i,j]=u_V[i,j]*a_V[i,j] 275 281 276 282 return … … 278 284 279 285 #-------------------------------------------------------- 280 #Boundaries - specific to the shallow_water_vel_domain286 #Boundaries - specific to the channel_domain 281 287 #-------------------------------------------------------- 282 288 class Reflective_boundary(Boundary): … … 298 304 #Handy shorthands 299 305 self.normals = domain.normals 300 self.area = domain.quantities['area'].vertex_values306 self.area = domain.quantities['area'].vertex_values 301 307 self.discharge = domain.quantities['discharge'].vertex_values 302 308 self.bed = domain.quantities['elevation'].vertex_values … … 306 312 self.stage = domain.quantities['stage'].vertex_values 307 313 308 from Numeric import zeros, Float 309 #self.conserved_quantities = zeros(3, Float) 310 self.evolved_quantities = zeros(7, Float) 314 self.evolved_quantities = numpy.zeros(7, numpy.float) 311 315 312 316 def __repr__(self): … … 320 324 321 325 q = self.evolved_quantities 322 q[0] = self.area[vol_id, edge_id]326 q[0] = self.area[vol_id, edge_id] 323 327 q[1] = -self.discharge[vol_id, edge_id] 324 q[2] = self.bed[vol_id, edge_id]325 q[3] = self.height[vol_id, edge_id]328 q[2] = self.bed[vol_id, edge_id] 329 q[3] = self.height[vol_id, edge_id] 326 330 q[4] = -self.velocity[vol_id, edge_id] 327 q[5] = self.width[vol_id,edge_id] 328 q[6] = self.stage[vol_id,edge_id] 329 330 #print "In Reflective q ",q 331 331 q[5] = self.width[vol_id,edge_id] 332 q[6] = self.stage[vol_id,edge_id] 332 333 333 334 return q … … 348 349 raise msg 349 350 350 from Numeric import array, Float 351 self.evolved_quantities=array(evolved_quantities).astype(Float) 351 assert len(evolved_quantities) == 7 352 353 self.evolved_quantities=numpy.array(evolved_quantities,numpy.float) 352 354 353 355 def __repr__(self): … … 370 372 371 373 372 Area = domain.quantities['area']374 Area = domain.quantities['area'] 373 375 Discharge = domain.quantities['discharge'] 374 376 Elevation = domain.quantities['elevation'] … … 402 404 def boundary_stress(domain): 403 405 404 405 406 from util import gradient 406 407 from Numeric import zeros, Float, array, sum 407 408 409 408 410 409 Area = domain.quantities['area'] … … 415 414 416 415 discharge_ud = Discharge.explicit_update 417 418 419 416 420 417 h = Height.vertex_values 421 418 b = Width.vertex_values … … 451 448 452 449 uh = domain.quantities['xmomentum'].centroid_values 453 #vh = domain.quantities['ymomentum'].centroid_values454 450 eta = domain.quantities['friction'].centroid_values 455 451 456 452 xmom_update = domain.quantities['xmomentum'].semi_implicit_update 457 #ymom_update = domain.quantities['ymomentum'].semi_implicit_update458 453 459 454 N = domain.number_of_elements … … 502 497 503 498 504 def check_forcefield(f):505 """Check that f is either506 1: a callable object f(t,x,y), where x and y are vectors507 and that it returns an array or a list of same length508 as x and y509 2: a scalar510 """511 512 from Numeric import ones, Float, array513 514 515 if callable(f):516 #N = 3517 N = 2518 #x = ones(3, Float)519 #y = ones(3, Float)520 x = ones(2, Float)521 #y = ones(2, Float)522 523 try:524 #q = f(1.0, x=x, y=y)525 q = f(1.0, x=x)526 except Exception, e:527 msg = 'Function %s could not be executed:\n%s' %(f, e)528 #FIXME: Reconsider this semantics529 raise msg530 531 try:532 q = array(q).astype(Float)533 except:534 msg = 'Return value from vector function %s could ' %f535 msg += 'not be converted into a Numeric array of floats.\n'536 msg += 'Specified function should return either list or array.'537 raise msg538 539 #Is this really what we want?540 msg = 'Return vector from function %s ' %f541 msg += 'must have same lenght as input vectors'542 assert len(q) == N, msg543 544 else:545 try:546 f = float(f)547 except:548 msg = 'Force field %s must be either a scalar' %f549 msg += ' or a vector function'550 raise msg551 return f552 499 553 500 def linearb(domain):
Note: See TracChangeset
for help on using the changeset viewer.