Changeset 7781 for anuga_work/development/pipeflow/sww_domain.py
 Timestamp:
 Jun 5, 2010, 7:01:28 PM (14 years ago)
 File:

 1 copied
Legend:
 Unmodified
 Added
 Removed

anuga_work/development/pipeflow/sww_domain.py
r7780 r7781 45 45 import numpy 46 46 47 from domain import * 48 Generic_Domain = Domain #Rename 47 from generic_domain import * 48 from sww_boundary_conditions import * 49 from sww_forcing_terms import * 49 50 50 51 #Shallow water domain … … 69 70 self.h0 = h0 70 71 71 #forcing terms not included in 1d domain ?WHy?72 #forcing terms 72 73 self.forcing_terms.append(gravity) 73 74 #self.forcing_terms.append(manning_friction) 74 #print "\nI have Removed forcing terms line 64 1dsw" 75 76 #Realtime visualisation 77 self.visualiser = None 78 self.visualise = False 79 self.visualise_color_stage = False 80 self.visualise_stage_range = 1.0 81 self.visualise_timer = True 82 self.visualise_range_z = None 75 83 76 84 77 #Stored output … … 131 124 132 125 133 def initialise_visualiser(self,scale_z=1.0,rect=None):134 #Realtime visualisation135 if self.visualiser is None:136 from realtime_visualisation_new import Visualiser137 self.visualiser = Visualiser(self,scale_z,rect)138 self.visualiser.setup['elevation']=True139 self.visualiser.updating['stage']=True140 self.visualise = True141 if self.visualise_color_stage == True:142 self.visualiser.coloring['stage'] = True143 self.visualiser.qcolor['stage'] = (0.0, 0.0, 0.8)144 print 'initialise visualiser'145 print self.visualiser.setup146 print self.visualiser.updating147 126 148 127 def check_integrity(self): … … 155 134 assert self.conserved_quantities[1] == 'xmomentum', msg 156 135 157 def extrapolate_second_order_sw(self): 158 #Call correct module function 159 #(either from this module or Cextension) 160 extrapolate_second_order_sw(self) 136 161 137 162 138 def compute_fluxes(self): 163 139 #Call correct module function 164 140 #(either from this module or Cextension) 165 compute_fluxes_C(self) 166 167 def compute_timestep(self): 168 #Call correct module function 169 compute_timestep(self) 141 142 import sys 143 144 145 timestep = numpy.float(sys.maxint) 146 147 stage = self.quantities['stage'] 148 xmom = self.quantities['xmomentum'] 149 bed = self.quantities['elevation'] 150 height = self.quantities['height'] 151 velocity = self.quantities['velocity'] 152 153 154 from sww_comp_flux_ext import compute_fluxes_ext 155 156 self.flux_timestep = compute_fluxes_ext(timestep,self,stage,xmom,bed,height,velocity) 157 170 158 171 159 def distribute_to_vertices_and_edges(self): … … 210 198 211 199 #=============== End of Shallow Water Domain =============================== 212 213 #Rotation of momentum vector214 def rotate(q, normal, direction = 1):215 """Rotate the momentum component q (q[1], q[2])216 from x,y coordinates to coordinates based on normal vector.217 218 If direction is negative the rotation is inverted.219 220 Input vector is preserved221 222 This function is specific to the shallow water wave equation223 """224 225 226 assert len(q) == 3,\227 'Vector of conserved quantities must have length 3'\228 'for 2D shallow water equation'229 230 try:231 l = len(normal)232 except:233 raise 'Normal vector must be an numpy array'234 235 assert l == 2, 'Normal vector must have 2 components'236 237 238 n1 = normal[0]239 n2 = normal[1]240 241 r = numpy.zeros(len(q), numpy.float) #Rotated quantities242 r[0] = q[0] #First quantity, height, is not rotated243 244 if direction == 1:245 n2 = n2246 247 248 r[1] = n1*q[1] + n2*q[2]249 r[2] = n2*q[1] + n1*q[2]250 251 return r252 253 254 def flux_function(normal, ql, qr, zl, zr):255 """Compute fluxes between volumes for the shallow water wave equation256 cast in terms of w = h+z using the 'central scheme' as described in257 258 Kurganov, Noelle, Petrova. 'Semidiscrete CentralUpwind Schemes For259 Hyperbolic Conservation Laws and HamiltonJacobi Equations'.260 Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707740.261 262 The implemented formula is given in equation (3.15) on page 714263 264 Conserved quantities w, uh, are stored as elements 0 and 1265 in the numpyal vectors ql an qr.266 267 Bed elevations zl and zr.268 """269 270 from config import g, epsilon, h0271 from math import sqrt272 273 #print 'ql',ql274 275 #Align momentums with xaxis276 #q_left = rotate(ql, normal, direction = 1)277 #q_right = rotate(qr, normal, direction = 1)278 q_left = ql279 q_left[1] = q_left[1]*normal280 q_right = qr281 q_right[1] = q_right[1]*normal282 283 #z = (zl+zr)/2 #Take average of field values284 z = 0.5*(zl+zr) #Take average of field values285 286 w_left = q_left[0] #w=h+z287 h_left = w_leftz288 uh_left = q_left[1]289 290 if h_left < epsilon:291 u_left = 0.0 #Could have been negative292 h_left = 0.0293 else:294 u_left = uh_left/(h_left + h0/h_left)295 296 297 uh_left = u_left*h_left298 299 300 w_right = q_right[0] #w=h+z301 h_right = w_rightz302 uh_right = q_right[1]303 304 if h_right < epsilon:305 u_right = 0.0 #Could have been negative306 h_right = 0.0307 else:308 u_right = uh_right/(h_right + h0/h_right)309 310 uh_right = u_right*h_right311 312 313 #We have got h and u at vertex, then the following is the calculation of fluxes!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!314 soundspeed_left = sqrt(g*h_left)315 soundspeed_right = sqrt(g*h_right)316 317 #Maximal wave speed318 s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)319 320 #Minimal wave speed321 s_min = min(u_leftsoundspeed_left, u_rightsoundspeed_right, 0)322 323 #Flux computation324 flux_left = numpy.array([u_left*h_left,325 u_left*uh_left + 0.5*g*h_left*h_left], numpy.float)326 flux_right = numpy.array([u_right*h_right,327 u_right*uh_right + 0.5*g*h_right*h_right], numpy.float)328 329 denom = s_maxs_min330 if denom == 0.0:331 edgeflux = array([0.0, 0.0])332 max_speed = 0.0333 else:334 edgeflux = (s_max*flux_left  s_min*flux_right)/denom335 edgeflux += s_max*s_min*(q_rightq_left)/denom336 337 edgeflux[1] = edgeflux[1]*normal338 339 max_speed = max(abs(s_max), abs(s_min))340 341 return edgeflux, max_speed342 #343 def compute_fluxes_python(domain):344 """345 Compute all fluxes and the timestep suitable for all volumes346 in domain.347 348 Compute total flux for each conserved quantity using "flux_function"349 350 Fluxes across each edge are scaled by edgelengths and summed up351 Resulting flux is then scaled by area and stored in352 explicit_update for each of the three conserved quantities353 stage, xmomentum and ymomentum354 355 The maximal allowable speed computed by the flux_function for each volume356 is converted to a timestep that must not be exceeded. The minimum of357 those is computed as the next overall timestep.358 359 Post conditions:360 domain.explicit_update is reset to computed flux values361 domain.timestep is set to the largest step satisfying all volumes.362 363 """364 365 import sys366 367 368 domain.distribute_to_vertices_and_edges()369 domain.update_boundary()370 371 N = domain.number_of_elements372 Stage = domain.quantities['stage']373 Xmom = domain.quantities['xmomentum']374 Bed = domain.quantities['elevation']375 376 stage = Stage.vertex_values377 xmom = Xmom.vertex_values378 bed = Bed.vertex_values379 380 stage_bdry = Stage.boundary_values381 xmom_bdry = Xmom.boundary_values382 383 384 385 flux = numpy.zeros(2, numpy.float) #Work array for summing up fluxes386 ql = numpy.zeros(2, numpy.float)387 qr = numpy.zeros(2, numpy.float)388 389 #Loop390 timestep = numpy.float(sys.maxint)391 enter = True392 for k in range(N):393 394 flux[:] = 0. #Reset work array395 #for i in range(3):396 for i in range(2):397 #Quantities inside volume facing neighbour i398 #ql[0] = stage[k, i]399 #ql[1] = xmom[k, i]400 ql = [stage[k, i], xmom[k, i]]401 zl = bed[k, i]402 403 #Quantities at neighbour on nearest face404 n = domain.neighbours[k,i]405 if n < 0:406 m = n1 #Convert negative flag to index407 qr[0] = stage_bdry[m]408 qr[1] = xmom_bdry[m]409 zr = zl #Extend bed elevation to boundary410 else:411 #m = domain.neighbour_edges[k,i]412 m = domain.neighbour_vertices[k,i]413 #qr = [stage[n, m], xmom[n, m], ymom[n, m]]414 qr[0] = stage[n, m]415 qr[1] = xmom[n, m]416 zr = bed[n, m]417 418 419 #Outward pointing normal vector420 normal = domain.normals[k, i]421 422 #Flux computation using provided function423 424 edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)425 426 #print 'edgeflux', edgeflux427 428 # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES429 # flux = edgefluxleft  edgefluxright430 flux = edgeflux #* domain.edgelengths[k,i]431 #Update optimal_timestep432 try:433 #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)434 timestep = min(timestep, domain.cfl*0.5*domain.areas[k]/max_speed)435 except ZeroDivisionError:436 pass437 438 #Normalise by area and store for when all conserved439 #quantities get updated440 flux /= domain.areas[k]441 442 Stage.explicit_update[k] = flux[0]443 Xmom.explicit_update[k] = flux[1]444 #Ymom.explicit_update[k] = flux[2]445 #print "flux cell",k,flux[0]446 447 domain.flux_timestep = timestep448 #print domain.quantities['stage'].centroid_values449 #450 def compute_timestep(domain):451 import sys452 453 454 N = domain.number_of_elements455 456 #Shortcuts457 Stage = domain.quantities['stage']458 Xmom = domain.quantities['xmomentum']459 Bed = domain.quantities['elevation']460 461 stage = Stage.vertex_values462 xmom = Xmom.vertex_values463 bed = Bed.vertex_values464 465 stage_bdry = Stage.boundary_values466 xmom_bdry = Xmom.boundary_values467 468 flux = zeros(2, numpy.float) #Work array for summing up fluxes469 ql = zeros(2, numpy.float)470 qr = zeros(2, numpy.float)471 472 #Loop473 timestep = numpy.float(sys.maxint)474 enter = True475 for k in range(N):476 477 flux[:] = 0. #Reset work array478 for i in range(2):479 #Quantities inside volume facing neighbour i480 ql = [stage[k, i], xmom[k, i]]481 zl = bed[k, i]482 483 #Quantities at neighbour on nearest face484 n = domain.neighbours[k,i]485 if n < 0:486 m = n1 #Convert negative flag to index487 qr[0] = stage_bdry[m]488 qr[1] = xmom_bdry[m]489 zr = zl #Extend bed elevation to boundary490 else:491 #m = domain.neighbour_edges[k,i]492 m = domain.neighbour_vertices[k,i]493 qr[0] = stage[n, m]494 qr[1] = xmom[n, m]495 zr = bed[n, m]496 497 498 #Outward pointing normal vector499 normal = domain.normals[k, i]500 501 edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)502 503 #Update optimal_timestep504 try:505 timestep = min(timestep, domain.cfl*0.5*domain.areas[k]/max_speed)506 except ZeroDivisionError:507 pass508 509 domain.timestep = timestep510 511 # Compute flux definition512 def compute_fluxes_C_long(domain):513 514 import sys515 516 517 timestep = numpy.float(sys.maxint)518 #print 'timestep=',timestep519 #print 'The type of timestep is',type(timestep)520 521 epsilon = domain.epsilon522 #print 'epsilon=',epsilon523 #print 'The type of epsilon is',type(epsilon)524 525 g = domain.g526 #print 'g=',g527 #print 'The type of g is',type(g)528 529 neighbours = domain.neighbours530 #print 'neighbours=',neighbours531 #print 'The type of neighbours is',type(neighbours)532 533 neighbour_vertices = domain.neighbour_vertices534 #print 'neighbour_vertices=',neighbour_vertices535 #print 'The type of neighbour_vertices is',type(neighbour_vertices)536 537 normals = domain.normals538 #print 'normals=',normals539 #print 'The type of normals is',type(normals)540 541 areas = domain.areas542 #print 'areas=',areas543 #print 'The type of areas is',type(areas)544 545 stage_edge_values = domain.quantities['stage'].vertex_values546 #print 'stage_edge_values=',stage_edge_values547 #print 'The type of stage_edge_values is',type(stage_edge_values)548 549 xmom_edge_values = domain.quantities['xmomentum'].vertex_values550 #print 'xmom_edge_values=',xmom_edge_values551 #print 'The type of xmom_edge_values is',type(xmom_edge_values)552 553 bed_edge_values = domain.quantities['elevation'].vertex_values554 #print 'bed_edge_values=',bed_edge_values555 #print 'The type of bed_edge_values is',type(bed_edge_values)556 557 stage_boundary_values = domain.quantities['stage'].boundary_values558 #print 'stage_boundary_values=',stage_boundary_values559 #print 'The type of stage_boundary_values is',type(stage_boundary_values)560 561 xmom_boundary_values = domain.quantities['xmomentum'].boundary_values562 #print 'xmom_boundary_values=',xmom_boundary_values563 #print 'The type of xmom_boundary_values is',type(xmom_boundary_values)564 565 stage_explicit_update = domain.quantities['stage'].explicit_update566 #print 'stage_explicit_update=',stage_explicit_update567 #print 'The type of stage_explicit_update is',type(stage_explicit_update)568 569 xmom_explicit_update = domain.quantities['xmomentum'].explicit_update570 #print 'xmom_explicit_update=',xmom_explicit_update571 #print 'The type of xmom_explicit_update is',type(xmom_explicit_update)572 573 number_of_elements = len(stage_edge_values)574 #print 'number_of_elements=',number_of_elements575 #print 'The type of number_of_elements is',type(number_of_elements)576 577 max_speed_array = domain.max_speed_array578 #print 'max_speed_array=',max_speed_array579 #print 'The type of max_speed_array is',type(max_speed_array)580 581 582 from comp_flux_ext import compute_fluxes_ext583 584 domain.flux_timestep = compute_fluxes_ext(timestep,585 epsilon,586 g,587 neighbours,588 neighbour_vertices,589 normals,590 areas,591 stage_edge_values,592 xmom_edge_values,593 bed_edge_values,594 stage_boundary_values,595 xmom_boundary_values,596 stage_explicit_update,597 xmom_explicit_update,598 number_of_elements,599 max_speed_array)600 601 602 # Compute flux definition603 def compute_fluxes_C(domain):604 import sys605 606 607 timestep = numpy.float(sys.maxint)608 609 stage = domain.quantities['stage']610 xmom = domain.quantities['xmomentum']611 bed = domain.quantities['elevation']612 height = domain.quantities['height']613 velocity = domain.quantities['velocity']614 615 616 from comp_flux_ext import compute_fluxes_ext_short617 618 domain.flux_timestep = compute_fluxes_ext_short(timestep,domain,stage,xmom,bed,height,velocity)619 620 621 622 # ###################################623 def compute_fluxes_C_wellbalanced(domain):624 #from numpy import zeros, numpy.float625 #import sys626 627 628 #timestep = numpy.float(sys.maxint)629 #epsilon = domain.epsilon630 #g = domain.g631 #neighbours = domain.neighbours632 #neighbour_vertices = domain.neighbour_vertices633 #normals = domain.normals634 #areas = domain.areas635 #stage_edge_values = domain.quantities['stage'].vertex_values636 #xmom_edge_values = domain.quantities['xmomentum'].vertex_values637 #bed_edge_values = domain.quantities['elevation'].vertex_values638 #stage_boundary_values = domain.quantities['stage'].boundary_values639 #xmom_boundary_values = domain.quantities['xmomentum'].boundary_values640 #stage_explicit_update = domain.quantities['stage'].explicit_update641 #xmom_explicit_update = domain.quantities['xmomentum'].explicit_update642 #number_of_elements = len(stage_edge_values)643 #max_speed_array = domain.max_speed_array644 645 import sys646 647 N = domain.number_of_elements648 timestep = numpy.float(sys.maxint)649 epsilon = domain.epsilon650 g = domain.g651 neighbours = domain.neighbours652 neighbour_vertices = domain.neighbour_vertices653 normals = domain.normals654 areas = domain.areas655 656 Stage = domain.quantities['stage']657 Xmom = domain.quantities['xmomentum']658 Bed = domain.quantities['elevation']659 660 stage_boundary_values = Stage.boundary_values661 xmom_boundary_values = Xmom.boundary_values662 stage_explicit_update = Stage.explicit_update663 xmom_explicit_update = Xmom.explicit_update664 max_speed_array = domain.max_speed_array665 666 domain.distribute_to_vertices_and_edges()667 domain.update_boundary()668 669 670 stage_V = Stage.vertex_values671 xmom_V = Xmom.vertex_values672 bed_V = Bed.vertex_values673 #h_V = Height.vertex_values674 #u_V = Velocity.vertex_values675 676 number_of_elements = len(stage_V)677 678 #flux = zeros(2, numpy.float) #Work array for summing up fluxes679 #ql = zeros(2, numpy.float)680 #qr = zeros(2, numpy.float)681 682 from comp_flux_ext_wellbalanced import compute_fluxes_ext_wellbalanced683 684 domain.flux_timestep = compute_fluxes_ext_wellbalanced(timestep,685 epsilon,686 g,687 neighbours,688 neighbour_vertices,689 normals,690 areas,691 stage_V,692 xmom_V,693 bed_V,694 stage_boundary_values,695 xmom_boundary_values,696 stage_explicit_update,697 xmom_explicit_update,698 number_of_elements,699 max_speed_array)700 701 # ###################################702 703 704 705 706 200 707 201 … … 1019 513 1020 514 1021 ############################################### 1022 #Boundaries  specific to the shallow water wave equation 1023 class Reflective_boundary(Boundary): 1024 """Reflective boundary returns same conserved quantities as 1025 those present in its neighbour volume but reflected. 1026 1027 This class is specific to the shallow water equation as it 1028 works with the momentum quantities assumed to be the second 1029 and third conserved quantities. 1030 """ 1031 1032 def __init__(self, domain = None): 1033 Boundary.__init__(self) 1034 1035 if domain is None: 1036 msg = 'Domain must be specified for reflective boundary' 1037 raise msg 1038 1039 #Handy shorthands 1040 self.normals = domain.normals 1041 self.stage = domain.quantities['stage'].vertex_values 1042 self.xmom = domain.quantities['xmomentum'].vertex_values 1043 self.bed = domain.quantities['elevation'].vertex_values 1044 self.height = domain.quantities['height'].vertex_values 1045 self.velocity = domain.quantities['velocity'].vertex_values 1046 1047 1048 self.quantities = numpy.zeros(5, numpy.float) 1049 1050 def __repr__(self): 1051 return 'Reflective_boundary' 1052 1053 1054 def evaluate(self, vol_id, edge_id): 1055 """Reflective boundaries reverses the outward momentum 1056 of the volume they serve. 1057 """ 1058 1059 q = self.quantities 1060 q[0] = self.stage[vol_id, edge_id] 1061 q[1] = self.xmom[vol_id, edge_id] 1062 q[2] = self.bed[vol_id, edge_id] 1063 q[3] = self.height[vol_id, edge_id] 1064 q[4] = self.velocity[vol_id, edge_id] 1065 1066 #normal = self.normals[vol_id,edge_id] 1067 1068 return q 1069 1070 class Dirichlet_boundary(Boundary): 1071 """Dirichlet boundary returns constant values for the 1072 conserved quantities 1073 """ 1074 1075 1076 def __init__(self, quantities=None): 1077 Boundary.__init__(self) 1078 1079 if quantities is None: 1080 msg = 'Must specify one value for each evolved quantity, w,uh,z,h,u' 1081 raise msg 1082 1083 self.quantities=numpy.array(quantities, numpy.float) 1084 1085 def __repr__(self): 1086 return 'Dirichlet boundary (%s)' %self.quantities 1087 1088 def evaluate(self, vol_id=None, edge_id=None): 1089 return self.quantities 1090 1091 1092 ######################### 1093 #Standard forcing terms: 1094 # 1095 def gravity(domain): 1096 """Apply gravitational pull in the presence of bed slope 1097 """ 1098 1099 from util import gradient 1100 1101 xmom = domain.quantities['xmomentum'].explicit_update 1102 stage = domain.quantities['stage'].explicit_update 1103 # ymom = domain.quantities['ymomentum'].explicit_update 1104 1105 Stage = domain.quantities['stage'] 1106 Elevation = domain.quantities['elevation'] 1107 #h = Stage.edge_values  Elevation.edge_values 1108 h = Stage.vertex_values  Elevation.vertex_values 1109 b = Elevation.vertex_values 1110 w = Stage.vertex_values 1111 1112 x = domain.get_vertex_coordinates() 1113 g = domain.g 1114 1115 for k in range(domain.number_of_elements): 1116 # avg_h = sum( h[k,:] )/3 1117 avg_h = sum( h[k,:] )/2 1118 1119 #Compute bed slope 1120 #x0, y0, x1, y1, x2, y2 = x[k,:] 1121 x0, x1 = x[k,:] 1122 #z0, z1, z2 = v[k,:] 1123 b0, b1 = b[k,:] 1124 1125 w0, w1 = w[k,:] 1126 wx = gradient(x0, x1, w0, w1) 1127 1128 #zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2) 1129 bx = gradient(x0, x1, b0, b1) 1130 1131 #Update momentum (explicit update is reset to source values) 1132 xmom[k] += g*bx*avg_h 1133 #xmom[k] = g*bx*avg_h 1134 #stage[k] = 0.0 1135 1136 1137 def manning_friction(domain): 1138 """Apply (Manning) friction to water momentum 1139 """ 1140 1141 from math import sqrt 1142 1143 w = domain.quantities['stage'].centroid_values 1144 z = domain.quantities['elevation'].centroid_values 1145 h = wz 1146 1147 uh = domain.quantities['xmomentum'].centroid_values 1148 #vh = domain.quantities['ymomentum'].centroid_values 1149 eta = domain.quantities['friction'].centroid_values 1150 1151 xmom_update = domain.quantities['xmomentum'].semi_implicit_update 1152 #ymom_update = domain.quantities['ymomentum'].semi_implicit_update 1153 1154 N = domain.number_of_elements 1155 eps = domain.minimum_allowed_height 1156 g = domain.g 1157 1158 for k in range(N): 1159 if eta[k] >= eps: 1160 if h[k] >= eps: 1161 #S = g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2)) 1162 S = g * eta[k]**2 * uh[k] 1163 S /= h[k]**(7.0/3) 1164 1165 #Update momentum 1166 xmom_update[k] += S*uh[k] 1167 #ymom_update[k] += S*vh[k] 1168 1169 def linear_friction(domain): 1170 """Apply linear friction to water momentum 1171 1172 Assumes quantity: 'linear_friction' to be present 1173 """ 1174 1175 from math import sqrt 1176 1177 w = domain.quantities['stage'].centroid_values 1178 z = domain.quantities['elevation'].centroid_values 1179 h = wz 1180 1181 uh = domain.quantities['xmomentum'].centroid_values 1182 # vh = domain.quantities['ymomentum'].centroid_values 1183 tau = domain.quantities['linear_friction'].centroid_values 1184 1185 xmom_update = domain.quantities['xmomentum'].semi_implicit_update 1186 # ymom_update = domain.quantities['ymomentum'].semi_implicit_update 1187 1188 N = domain.number_of_elements 1189 eps = domain.minimum_allowed_height 1190 g = domain.g #Not necessary? Why was this added? 1191 1192 for k in range(N): 1193 if tau[k] >= eps: 1194 if h[k] >= eps: 1195 S = tau[k]/h[k] 1196 1197 #Update momentum 1198 xmom_update[k] += S*uh[k] 1199 # ymom_update[k] += S*vh[k] 1200 1201 1202 1203 def check_forcefield(f): 1204 """Check that f is either 1205 1: a callable object f(t,x,y), where x and y are vectors 1206 and that it returns an array or a list of same length 1207 as x and y 1208 2: a scalar 1209 """ 1210 1211 1212 if callable(f): 1213 #N = 3 1214 N = 2 1215 #x = ones(3, numpy.float) 1216 #y = ones(3, numpy.float) 1217 x = ones(2, numpy.float) 1218 #y = ones(2, numpy.float) 1219 1220 try: 1221 #q = f(1.0, x=x, y=y) 1222 q = f(1.0, x=x) 1223 except Exception, e: 1224 msg = 'Function %s could not be executed:\n%s' %(f, e) 1225 #FIXME: Reconsider this semantics 1226 raise msg 1227 1228 try: 1229 q = numpy.array(q, numpy.float) 1230 except: 1231 msg = 'Return value from vector function %s could ' %f 1232 msg += 'not be converted into a numpy array of numpy.floats.\n' 1233 msg += 'Specified function should return either list or array.' 1234 raise msg 1235 1236 #Is this really what we want? 1237 msg = 'Return vector from function %s ' %f 1238 msg += 'must have same lenght as input vectors' 1239 assert len(q) == N, msg 1240 1241 else: 1242 try: 1243 f = numpy.float(f) 1244 except: 1245 msg = 'Force field %s must be either a scalar' %f 1246 msg += ' or a vector function' 1247 raise msg 1248 return f 1249 1250 class Wind_stress: 1251 """Apply wind stress to water momentum in terms of 1252 wind speed [m/s] and wind direction [degrees] 1253 """ 1254 1255 def __init__(self, *args, **kwargs): 1256 """Initialise windfield from wind speed s [m/s] 1257 and wind direction phi [degrees] 1258 1259 Inputs v and phi can be either scalars or Python functions, e.g. 1260 1261 W = Wind_stress(10, 178) 1262 1263 #FIXME  'normal' degrees are assumed for now, i.e. the 1264 vector (1,0) has zero degrees. 1265 We may need to convert from 'compass' degrees later on and also 1266 map from True north to grid north. 1267 1268 Arguments can also be Python functions of t,x,y as in 1269 1270 def speed(t,x,y): 1271 ... 1272 return s 1273 1274 def angle(t,x,y): 1275 ... 1276 return phi 1277 1278 where x and y are vectors. 1279 1280 and then pass the functions in 1281 1282 W = Wind_stress(speed, angle) 1283 1284 The instantiated object W can be appended to the list of 1285 forcing_terms as in 1286 1287 Alternatively, one vector valued function for (speed, angle) 1288 can be applied, providing both quantities simultaneously. 1289 As in 1290 W = Wind_stress(F), where returns (speed, angle) for each t. 1291 1292 domain.forcing_terms.append(W) 1293 """ 1294 1295 from config import rho_a, rho_w, eta_w 1296 1297 if len(args) == 2: 1298 s = args[0] 1299 phi = args[1] 1300 elif len(args) == 1: 1301 #Assume vector function returning (s, phi)(t,x,y) 1302 vector_function = args[0] 1303 #s = lambda t,x,y: vector_function(t,x=x,y=y)[0] 1304 #phi = lambda t,x,y: vector_function(t,x=x,y=y)[1] 1305 s = lambda t,x: vector_function(t,x=x)[0] 1306 phi = lambda t,x: vector_function(t,x=x)[1] 1307 else: 1308 #Assume info is in 2 keyword arguments 1309 1310 if len(kwargs) == 2: 1311 s = kwargs['s'] 1312 phi = kwargs['phi'] 1313 else: 1314 raise 'Assumes two keyword arguments: s=..., phi=....' 1315 1316 print 'phi', phi 1317 self.speed = check_forcefield(s) 1318 self.phi = check_forcefield(phi) 1319 1320 self.const = eta_w*rho_a/rho_w 1321 1322 1323 def __call__(self, domain): 1324 """Evaluate windfield based on values found in domain 1325 """ 1326 1327 from math import pi, cos, sin, sqrt 1328 1329 xmom_update = domain.quantities['xmomentum'].explicit_update 1330 #ymom_update = domain.quantities['ymomentum'].explicit_update 1331 1332 N = domain.number_of_elements 1333 t = domain.time 1334 1335 if callable(self.speed): 1336 xc = domain.get_centroid_coordinates() 1337 #s_vec = self.speed(t, xc[:,0], xc[:,1]) 1338 s_vec = self.speed(t, xc) 1339 else: 1340 #Assume s is a scalar 1341 1342 try: 1343 s_vec = self.speed * ones(N, numpy.float) 1344 except: 1345 msg = 'Speed must be either callable or a scalar: %s' %self.s 1346 raise msg 1347 1348 1349 if callable(self.phi): 1350 xc = domain.get_centroid_coordinates() 1351 #phi_vec = self.phi(t, xc[:,0], xc[:,1]) 1352 phi_vec = self.phi(t, xc) 1353 else: 1354 #Assume phi is a scalar 1355 1356 try: 1357 phi_vec = self.phi * ones(N, numpy.float) 1358 except: 1359 msg = 'Angle must be either callable or a scalar: %s' %self.phi 1360 raise msg 1361 1362 #assign_windfield_values(xmom_update, ymom_update, 1363 # s_vec, phi_vec, self.const) 1364 assign_windfield_values(xmom_update, s_vec, phi_vec, self.const) 1365 1366 1367 #def assign_windfield_values(xmom_update, ymom_update, 1368 # s_vec, phi_vec, const): 1369 def assign_windfield_values(xmom_update, s_vec, phi_vec, const): 1370 """Python version of assigning wind field to update vectors. 1371 A c version also exists (for speed) 1372 """ 1373 from math import pi, cos, sin, sqrt 1374 1375 N = len(s_vec) 1376 for k in range(N): 1377 s = s_vec[k] 1378 phi = phi_vec[k] 1379 1380 #Convert to radians 1381 phi = phi*pi/180 1382 1383 #Compute velocity vector (u, v) 1384 u = s*cos(phi) 1385 v = s*sin(phi) 1386 1387 #Compute wind stress 1388 #S = const * sqrt(u**2 + v**2) 1389 S = const * u 1390 xmom_update[k] += S*u 1391 #ymom_update[k] += S*v 515
Note: See TracChangeset
for help on using the changeset viewer.