"""Class Domain - 1D interval domains for finite-volume computations of the shallow water wave equation. This module contains a specialisation of class Domain from module domain.py consisting of methods specific to the Shallow Water Wave Equation U_t + E_x = S where U = [w, uh] E = [uh, u^2h + gh^2/2] S represents source terms forcing the system (e.g. gravity, friction, wind stress, ...) and _t, _x, _y denote the derivative with respect to t, x and y respectiely. The quantities are symbol variable name explanation x x horizontal distance from origin [m] z elevation elevation of bed on which flow is modelled [m] h height water height above z [m] w stage absolute water level, w = z+h [m] u speed in the x direction [m/s] uh xmomentum momentum in the x direction [m^2/s] eta mannings friction coefficient [to appear] nu wind stress coefficient [to appear] The conserved quantities are w, uh For details see e.g. Christopher Zoppou and Stephen Roberts, Catastrophic Collapse of Water Supply Reservoirs in Urban Areas, Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999 John Jakeman, Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou Geoscience Australia, 2006 """ #from domain import * #from domain_order2 import * from domain_t2 import * from comp_flux_ext import * Generic_Domain = Domain #Rename #Shallow water domain class Domain(Generic_Domain): def __init__(self, coordinates, boundary = None, tagged_elements = None, geo_reference = None): conserved_quantities = ['stage', 'xmomentum'] other_quantities = ['elevation', 'friction'] Generic_Domain.__init__(self, coordinates, boundary, conserved_quantities, other_quantities, tagged_elements, geo_reference) from config import minimum_allowed_height, g self.minimum_allowed_height = minimum_allowed_height self.g = g #forcing terms not included in 1d domain ?WHy? self.forcing_terms.append(gravity) #self.forcing_terms.append(manning_friction) #print "\nI have Removed forcing terms line 64 1dsw" #Realtime visualisation self.visualiser = None self.visualise = False self.visualise_color_stage = False self.visualise_stage_range = 1.0 self.visualise_timer = True self.visualise_range_z = None #Stored output self.store = True self.format = 'sww' self.smooth = True #Evolve parametrs self.cfl = 1.0 #Reduction operation for get_vertex_values from util import mean self.reduction = mean #self.reduction = min #Looks better near steep slopes self.quantities_to_be_stored = ['stage','xmomentum'] def set_quantities_to_be_stored(self, q): """Specify which quantities will be stored in the sww file. q must be either: - the name of a quantity - a list of quantity names - None In the two first cases, the named quantities will be stored at each yieldstep (This is in addition to the quantities elevation and friction) If q is None, storage will be switched off altogether. """ if q is None: self.quantities_to_be_stored = [] self.store = False return if isinstance(q, basestring): q = [q] # Turn argument into a list #Check correcness for quantity_name in q: msg = 'Quantity %s is not a valid conserved quantity' %quantity_name assert quantity_name in self.conserved_quantities, msg self.quantities_to_be_stored = q def initialise_visualiser(self,scale_z=1.0,rect=None): #Realtime visualisation if self.visualiser is None: from realtime_visualisation_new import Visualiser self.visualiser = Visualiser(self,scale_z,rect) self.visualiser.setup['elevation']=True self.visualiser.updating['stage']=True self.visualise = True if self.visualise_color_stage == True: self.visualiser.coloring['stage'] = True self.visualiser.qcolor['stage'] = (0.0, 0.0, 0.8) print 'initialise visualiser' print self.visualiser.setup print self.visualiser.updating def check_integrity(self): Generic_Domain.check_integrity(self) #Check that we are solving the shallow water wave equation msg = 'First conserved quantity must be "stage"' assert self.conserved_quantities[0] == 'stage', msg msg = 'Second conserved quantity must be "xmomentum"' assert self.conserved_quantities[1] == 'xmomentum', msg def extrapolate_second_order_sw(self): #Call correct module function #(either from this module or C-extension) extrapolate_second_order_sw(self) def compute_fluxes(self): #Call correct module function #(either from this module or C-extension) compute_fluxes(self) def compute_timestep(self): #Call correct module function compute_timestep(self) def distribute_to_vertices_and_edges(self): #Call correct module function #(either from this module or C-extension) distribute_to_vertices_and_edges(self) def evolve(self, yieldstep = None, finaltime = None, skip_initial_step = False): """Specialisation of basic evolve method from parent class """ #Call check integrity here rather than from user scripts #self.check_integrity() #msg = 'Parameter beta_h must be in the interval [0, 1)' #assert 0 <= self.beta_h < 1.0, msg #msg = 'Parameter beta_w must be in the interval [0, 1)' #assert 0 <= self.beta_w < 1.0, msg #Initial update of vertex and edge values before any storage #and or visualisation self.distribute_to_vertices_and_edges() #Call basic machinery from parent class for t in Generic_Domain.evolve(self, yieldstep, finaltime, skip_initial_step): yield(t) def initialise_storage(self): """Create and initialise self.writer object for storing data. Also, save x and bed elevation """ import data_manager #Initialise writer self.writer = data_manager.get_dataobject(self, mode = 'w') #Store vertices and connectivity self.writer.store_connectivity() def store_timestep(self, name): """Store named quantity and time. Precondition: self.write has been initialised """ self.writer.store_timestep(name) #=============== End of Shallow Water Domain =============================== #Rotation of momentum vector def rotate(q, normal, direction = 1): """Rotate the momentum component q (q[1], q[2]) from x,y coordinates to coordinates based on normal vector. If direction is negative the rotation is inverted. Input vector is preserved This function is specific to the shallow water wave equation """ from Numeric import zeros, Float assert len(q) == 3,\ 'Vector of conserved quantities must have length 3'\ 'for 2D shallow water equation' try: l = len(normal) except: raise 'Normal vector must be an Numeric array' assert l == 2, 'Normal vector must have 2 components' n1 = normal[0] n2 = normal[1] r = zeros(len(q), Float) #Rotated quantities r[0] = q[0] #First quantity, height, is not rotated if direction == -1: n2 = -n2 r[1] = n1*q[1] + n2*q[2] r[2] = -n2*q[1] + n1*q[2] return r """++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++""" # Compute flux definition def compute_fluxes(domain): from Numeric import zeros, Float import sys timestep = float(sys.maxint) print 'timestep=',timestep print 'The type of timestep is',type(timestep) epsilon = domain.epsilon print 'epsilon=',epsilon print 'The type of epsilon is',type(epsilon) g = domain.g print 'g=',g print 'The type of g is',type(g) neighbours = domain.neighbours print 'neighbours=',neighbours print 'The type of neighbours is',type(neighbours) neighbour_vertices = domain.neighbour_vertices print 'neighbour_vertices=',neighbour_vertices print 'The type of neighbour_vertices is',type(neighbour_vertices) normals = domain.normals print 'normals=',normals print 'The type of normals is',type(normals) areas = domain.areas print 'areas=',areas print 'The type of areas is',type(areas) stage_edge_values = domain.quantities['stage'].vertex_values print 'stage_edge_values=',stage_edge_values print 'The type of stage_edge_values is',type(stage_edge_values) xmom_edge_values = domain.quantities['xmomentum'].vertex_values print 'xmom_edge_values=',xmom_edge_values print 'The type of xmom_edge_values is',type(xmom_edge_values) bed_edge_values = domain.quantities['elevation'].vertex_values print 'bed_edge_values=',bed_edge_values print 'The type of bed_edge_values is',type(bed_edge_values) stage_boundary_values = domain.quantities['stage'].boundary_values print 'stage_boundary_values=',stage_boundary_values print 'The type of stage_boundary_values is',type(stage_boundary_values) xmom_boundary_values = domain.quantities['xmomentum'].boundary_values print 'xmom_boundary_values=',xmom_boundary_values print 'The type of xmom_boundary_values is',type(xmom_boundary_values) stage_explicit_update = domain.quantities['stage'].explicit_update print 'stage_explicit_update=',stage_explicit_update print 'The type of stage_explicit_update is',type(stage_explicit_update) xmom_explicit_update = domain.quantities['xmomentum'].explicit_update print 'xmom_explicit_update=',xmom_explicit_update print 'The type of xmom_explicit_update is',type(xmom_explicit_update) number_of_elements = len(stage_edge_values) print 'number_of_elements=',number_of_elements print 'The type of number_of_elements is',type(number_of_elements) max_speed_array = zeros(number_of_elements, Float) print 'max_speed_array=',max_speed_array print 'The type of max_speed_array is',type(max_speed_array) return compute_fluxes_ext(timestep, epsilon, g, neighbours, neighbour_vertices, normals, areas, stage_edge_values, xmom_edge_values, bed_edge_values, stage_boundary_values, xmom_boundary_values, number_of_elements, max_speed_array) #################################### # Module functions for gradient limiting (distribute_to_vertices_and_edges) def distribute_to_vertices_and_edges(domain): """Distribution from centroids to vertices specific to the shallow water wave equation. It will ensure that h (w-z) is always non-negative even in the presence of steep bed-slopes by taking a weighted average between shallow and deep cases. In addition, all conserved quantities get distributed as per either a constant (order==1) or a piecewise linear function (order==2). FIXME: more explanation about removal of artificial variability etc Precondition: All quantities defined at centroids and bed elevation defined at vertices. Postcondition Conserved quantities defined at vertices """ #from config import optimised_gradient_limiter #Remove very thin layers of water protect_against_infinitesimal_and_negative_heights(domain) #Extrapolate all conserved quantities #if optimised_gradient_limiter: # #MH090605 if second order, # #perform the extrapolation and limiting on # #all of the conserved quantities # if (domain.order == 1): # for name in domain.conserved_quantities: # Q = domain.quantities[name] # Q.extrapolate_first_order() # elif domain.order == 2: # domain.extrapolate_second_order_sw() # else: # raise 'Unknown order' #else: #old code: for name in domain.conserved_quantities: Q = domain.quantities[name] if domain.order == 1: Q.extrapolate_first_order() elif domain.order == 2: #print "add extrapolate second order to shallow water" #if name != 'height': Q.extrapolate_second_order() #Q.limit() else: raise 'Unknown order' #Take bed elevation into account when water heights are small #balance_deep_and_shallow(domain) #protect_against_infinitesimal_and_negative_heights(domain) #Compute edge values by interpolation #for name in domain.conserved_quantities: # Q = domain.quantities[name] # Q.interpolate_from_vertices_to_edges() def protect_against_infinitesimal_and_negative_heights(domain): """Protect against infinitesimal heights and associated high velocities """ #Shortcuts wc = domain.quantities['stage'].centroid_values zc = domain.quantities['elevation'].centroid_values xmomc = domain.quantities['xmomentum'].centroid_values # ymomc = domain.quantities['ymomentum'].centroid_values hc = wc - zc #Water depths at centroids zv = domain.quantities['elevation'].vertex_values wv = domain.quantities['stage'].vertex_values hv = wv-zv xmomv = domain.quantities['xmomentum'].vertex_values #remove the above two lines and corresponding code below #Update for k in range(domain.number_of_elements): if hc[k] < domain.minimum_allowed_height: #Control stage if hc[k] < domain.epsilon: wc[k] = zc[k] # Contain 'lost mass' error wv[k,0] = zv[k,0] wv[k,1] = zv[k,1] xmomc[k] = 0.0 #N = domain.number_of_elements #if (k == 0) | (k==N-1): # wc[k] = zc[k] # Contain 'lost mass' error # wv[k,0] = zv[k,0] # wv[k,1] = zv[k,1] def h_limiter(domain): """Limit slopes for each volume to eliminate artificial variance introduced by e.g. second order extrapolator limit on h = w-z This limiter depends on two quantities (w,z) so it resides within this module rather than within quantity.py """ from Numeric import zeros, Float N = domain.number_of_elements beta_h = domain.beta_h #Shortcuts wc = domain.quantities['stage'].centroid_values zc = domain.quantities['elevation'].centroid_values hc = wc - zc wv = domain.quantities['stage'].vertex_values zv = domain.quantities['elevation'].vertex_values hv = wv-zv hvbar = zeros(hv.shape, Float) #h-limited values #Find min and max of this and neighbour's centroid values hmax = zeros(hc.shape, Float) hmin = zeros(hc.shape, Float) for k in range(N): hmax[k] = hmin[k] = hc[k] #for i in range(3): for i in range(2): n = domain.neighbours[k,i] if n >= 0: hn = hc[n] #Neighbour's centroid value hmin[k] = min(hmin[k], hn) hmax[k] = max(hmax[k], hn) #Diffences between centroids and maxima/minima dhmax = hmax - hc dhmin = hmin - hc #Deltas between vertex and centroid values dh = zeros(hv.shape, Float) #for i in range(3): for i in range(2): dh[:,i] = hv[:,i] - hc #Phi limiter for k in range(N): #Find the gradient limiter (phi) across vertices phi = 1.0 #for i in range(3): for i in range(2): r = 1.0 if (dh[k,i] > 0): r = dhmax[k]/dh[k,i] if (dh[k,i] < 0): r = dhmin[k]/dh[k,i] phi = min( min(r*beta_h, 1), phi ) #Then update using phi limiter #for i in range(3): for i in range(2): hvbar[k,i] = hc[k] + phi*dh[k,i] return hvbar def balance_deep_and_shallow(domain): """Compute linear combination between stage as computed by gradient-limiters limiting using w, and stage computed by gradient-limiters limiting using h (h-limiter). The former takes precedence when heights are large compared to the bed slope while the latter takes precedence when heights are relatively small. Anything in between is computed as a balanced linear combination in order to avoid numerical disturbances which would otherwise appear as a result of hard switching between modes. The h-limiter is always applied irrespective of the order. """ #Shortcuts wc = domain.quantities['stage'].centroid_values zc = domain.quantities['elevation'].centroid_values hc = wc - zc wv = domain.quantities['stage'].vertex_values zv = domain.quantities['elevation'].vertex_values hv = wv-zv #Limit h hvbar = h_limiter(domain) for k in range(domain.number_of_elements): #Compute maximal variation in bed elevation # This quantitiy is # dz = max_i abs(z_i - z_c) # and it is independent of dimension # In the 1d case zc = (z0+z1)/2 # In the 2d case zc = (z0+z1+z2)/3 dz = max(abs(zv[k,0]-zc[k]), abs(zv[k,1]-zc[k]))#, # abs(zv[k,2]-zc[k])) hmin = min( hv[k,:] ) #Create alpha in [0,1], where alpha==0 means using the h-limited #stage and alpha==1 means using the w-limited stage as #computed by the gradient limiter (both 1st or 2nd order) #If hmin > dz/2 then alpha = 1 and the bed will have no effect #If hmin < 0 then alpha = 0 reverting to constant height above bed. if dz > 0.0: alpha = max( min( 2*hmin/dz, 1.0), 0.0 ) else: #Flat bed alpha = 1.0 alpha = 0.0 #Let # # wvi be the w-limited stage (wvi = zvi + hvi) # wvi- be the h-limited state (wvi- = zvi + hvi-) # # #where i=0,1,2 denotes the vertex ids # #Weighted balance between w-limited and h-limited stage is # # wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi) # #It follows that the updated wvi is # wvi := zvi + (1-alpha)*hvi- + alpha*hvi # # Momentum is balanced between constant and limited #for i in range(3): # wv[k,i] = zv[k,i] + hvbar[k,i] #return if alpha < 1: #for i in range(3): for i in range(2): wv[k,i] = zv[k,i] + (1.0-alpha)*hvbar[k,i] + alpha*hv[k,i] #Momentums at centroids xmomc = domain.quantities['xmomentum'].centroid_values # ymomc = domain.quantities['ymomentum'].centroid_values #Momentums at vertices xmomv = domain.quantities['xmomentum'].vertex_values # ymomv = domain.quantities['ymomentum'].vertex_values # Update momentum as a linear combination of # xmomc and ymomc (shallow) and momentum # from extrapolator xmomv and ymomv (deep). xmomv[k,:] = (1.0-alpha)*xmomc[k] + alpha*xmomv[k,:] # ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:] ############################################### #Boundaries - specific to the shallow water wave equation class Reflective_boundary(Boundary): """Reflective boundary returns same conserved quantities as those present in its neighbour volume but reflected. This class is specific to the shallow water equation as it works with the momentum quantities assumed to be the second and third conserved quantities. """ def __init__(self, domain = None): Boundary.__init__(self) if domain is None: msg = 'Domain must be specified for reflective boundary' raise msg #Handy shorthands #self.stage = domain.quantities['stage'].edge_values #self.xmom = domain.quantities['xmomentum'].edge_values #self.ymom = domain.quantities['ymomentum'].edge_values self.normals = domain.normals self.stage = domain.quantities['stage'].vertex_values self.xmom = domain.quantities['xmomentum'].vertex_values from Numeric import zeros, Float #self.conserved_quantities = zeros(3, Float) self.conserved_quantities = zeros(2, Float) def __repr__(self): return 'Reflective_boundary' def evaluate(self, vol_id, edge_id): """Reflective boundaries reverses the outward momentum of the volume they serve. """ q = self.conserved_quantities q[0] = self.stage[vol_id, edge_id] q[1] = self.xmom[vol_id, edge_id] #q[2] = self.ymom[vol_id, edge_id] #normal = self.normals[vol_id, 2*edge_id:2*edge_id+2] #normal = self.normals[vol_id, 2*edge_id:2*edge_id+1] normal = self.normals[vol_id,edge_id] #r = rotate(q, normal, direction = 1) #r[1] = -r[1] #q = rotate(r, normal, direction = -1) r = q r[1] = normal*r[1] r[1] = -r[1] r[1] = normal*r[1] q = r #For start interval there is no outward momentum so do not need to #reverse direction in this case return q class Dirichlet_boundary(Boundary): """Dirichlet boundary returns constant values for the conserved quantities """ def __init__(self, conserved_quantities=None): Boundary.__init__(self) if conserved_quantities is None: msg = 'Must specify one value for each conserved quantity' raise msg from Numeric import array, Float self.conserved_quantities=array(conserved_quantities).astype(Float) def __repr__(self): return 'Dirichlet boundary (%s)' %self.conserved_quantities def evaluate(self, vol_id=None, edge_id=None): return self.conserved_quantities ######################### #Standard forcing terms: # def gravity(domain): """Apply gravitational pull in the presence of bed slope """ from util import gradient from Numeric import zeros, Float, array, sum xmom = domain.quantities['xmomentum'].explicit_update stage = domain.quantities['stage'].explicit_update # ymom = domain.quantities['ymomentum'].explicit_update Stage = domain.quantities['stage'] Elevation = domain.quantities['elevation'] #h = Stage.edge_values - Elevation.edge_values h = Stage.vertex_values - Elevation.vertex_values b = Elevation.vertex_values w = Stage.vertex_values x = domain.get_vertex_coordinates() g = domain.g for k in range(domain.number_of_elements): # avg_h = sum( h[k,:] )/3 avg_h = sum( h[k,:] )/2 #Compute bed slope #x0, y0, x1, y1, x2, y2 = x[k,:] x0, x1 = x[k,:] #z0, z1, z2 = v[k,:] b0, b1 = b[k,:] w0, w1 = w[k,:] wx = gradient(x0, x1, w0, w1) #zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2) bx = gradient(x0, x1, b0, b1) #Update momentum (explicit update is reset to source values) if domain.split == False: xmom[k] += -g*bx*avg_h #xmom[k] = -g*bx*avg_h #stage[k] = 0.0 elif domain.split == True: xmom[k] += -g*wx*avg_h #xmom[k] = -g*wx*avg_h #ymom[k] += -g*zy*avg_h def manning_friction(domain): """Apply (Manning) friction to water momentum """ from math import sqrt w = domain.quantities['stage'].centroid_values z = domain.quantities['elevation'].centroid_values h = w-z uh = domain.quantities['xmomentum'].centroid_values #vh = domain.quantities['ymomentum'].centroid_values eta = domain.quantities['friction'].centroid_values xmom_update = domain.quantities['xmomentum'].semi_implicit_update #ymom_update = domain.quantities['ymomentum'].semi_implicit_update N = domain.number_of_elements eps = domain.minimum_allowed_height g = domain.g for k in range(N): if eta[k] >= eps: if h[k] >= eps: #S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2)) S = -g * eta[k]**2 * uh[k] S /= h[k]**(7.0/3) #Update momentum xmom_update[k] += S*uh[k] #ymom_update[k] += S*vh[k] def linear_friction(domain): """Apply linear friction to water momentum Assumes quantity: 'linear_friction' to be present """ from math import sqrt w = domain.quantities['stage'].centroid_values z = domain.quantities['elevation'].centroid_values h = w-z uh = domain.quantities['xmomentum'].centroid_values # vh = domain.quantities['ymomentum'].centroid_values tau = domain.quantities['linear_friction'].centroid_values xmom_update = domain.quantities['xmomentum'].semi_implicit_update # ymom_update = domain.quantities['ymomentum'].semi_implicit_update N = domain.number_of_elements eps = domain.minimum_allowed_height g = domain.g #Not necessary? Why was this added? for k in range(N): if tau[k] >= eps: if h[k] >= eps: S = -tau[k]/h[k] #Update momentum xmom_update[k] += S*uh[k] # ymom_update[k] += S*vh[k] def check_forcefield(f): """Check that f is either 1: a callable object f(t,x,y), where x and y are vectors and that it returns an array or a list of same length as x and y 2: a scalar """ from Numeric import ones, Float, array if callable(f): #N = 3 N = 2 #x = ones(3, Float) #y = ones(3, Float) x = ones(2, Float) #y = ones(2, Float) try: #q = f(1.0, x=x, y=y) q = f(1.0, x=x) except Exception, e: msg = 'Function %s could not be executed:\n%s' %(f, e) #FIXME: Reconsider this semantics raise msg try: q = array(q).astype(Float) except: msg = 'Return value from vector function %s could ' %f msg += 'not be converted into a Numeric array of floats.\n' msg += 'Specified function should return either list or array.' raise msg #Is this really what we want? msg = 'Return vector from function %s ' %f msg += 'must have same lenght as input vectors' assert len(q) == N, msg else: try: f = float(f) except: msg = 'Force field %s must be either a scalar' %f msg += ' or a vector function' raise msg return f class Wind_stress: """Apply wind stress to water momentum in terms of wind speed [m/s] and wind direction [degrees] """ def __init__(self, *args, **kwargs): """Initialise windfield from wind speed s [m/s] and wind direction phi [degrees] Inputs v and phi can be either scalars or Python functions, e.g. W = Wind_stress(10, 178) #FIXME - 'normal' degrees are assumed for now, i.e. the vector (1,0) has zero degrees. We may need to convert from 'compass' degrees later on and also map from True north to grid north. Arguments can also be Python functions of t,x,y as in def speed(t,x,y): ... return s def angle(t,x,y): ... return phi where x and y are vectors. and then pass the functions in W = Wind_stress(speed, angle) The instantiated object W can be appended to the list of forcing_terms as in Alternatively, one vector valued function for (speed, angle) can be applied, providing both quantities simultaneously. As in W = Wind_stress(F), where returns (speed, angle) for each t. domain.forcing_terms.append(W) """ from config import rho_a, rho_w, eta_w from Numeric import array, Float if len(args) == 2: s = args[0] phi = args[1] elif len(args) == 1: #Assume vector function returning (s, phi)(t,x,y) vector_function = args[0] #s = lambda t,x,y: vector_function(t,x=x,y=y)[0] #phi = lambda t,x,y: vector_function(t,x=x,y=y)[1] s = lambda t,x: vector_function(t,x=x)[0] phi = lambda t,x: vector_function(t,x=x)[1] else: #Assume info is in 2 keyword arguments if len(kwargs) == 2: s = kwargs['s'] phi = kwargs['phi'] else: raise 'Assumes two keyword arguments: s=..., phi=....' print 'phi', phi self.speed = check_forcefield(s) self.phi = check_forcefield(phi) self.const = eta_w*rho_a/rho_w def __call__(self, domain): """Evaluate windfield based on values found in domain """ from math import pi, cos, sin, sqrt from Numeric import Float, ones, ArrayType xmom_update = domain.quantities['xmomentum'].explicit_update #ymom_update = domain.quantities['ymomentum'].explicit_update N = domain.number_of_elements t = domain.time if callable(self.speed): xc = domain.get_centroid_coordinates() #s_vec = self.speed(t, xc[:,0], xc[:,1]) s_vec = self.speed(t, xc) else: #Assume s is a scalar try: s_vec = self.speed * ones(N, Float) except: msg = 'Speed must be either callable or a scalar: %s' %self.s raise msg if callable(self.phi): xc = domain.get_centroid_coordinates() #phi_vec = self.phi(t, xc[:,0], xc[:,1]) phi_vec = self.phi(t, xc) else: #Assume phi is a scalar try: phi_vec = self.phi * ones(N, Float) except: msg = 'Angle must be either callable or a scalar: %s' %self.phi raise msg #assign_windfield_values(xmom_update, ymom_update, # s_vec, phi_vec, self.const) assign_windfield_values(xmom_update, s_vec, phi_vec, self.const) #def assign_windfield_values(xmom_update, ymom_update, # s_vec, phi_vec, const): def assign_windfield_values(xmom_update, s_vec, phi_vec, const): """Python version of assigning wind field to update vectors. A c version also exists (for speed) """ from math import pi, cos, sin, sqrt N = len(s_vec) for k in range(N): s = s_vec[k] phi = phi_vec[k] #Convert to radians phi = phi*pi/180 #Compute velocity vector (u, v) u = s*cos(phi) v = s*sin(phi) #Compute wind stress #S = const * sqrt(u**2 + v**2) S = const * u xmom_update[k] += S*u #ymom_update[k] += S*v