"""Class Domain - 2D triangular 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 FIXME: Write equations here! Conserved quantities are w (water level or stage), uh (x momentum) and vh (y momentum). Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou Geoscience Australia, 2004 """ from domain import * Generic_domain = Domain #Rename class Domain(Generic_domain): def __init__(self, coordinates, vertices, boundary = None): conserved_quantities = ['level', 'xmomentum', 'ymomentum'] other_quantities = ['elevation', 'friction'] Generic_domain.__init__(self, coordinates, vertices, boundary, conserved_quantities, other_quantities) from config import minimum_allowed_height self.minimum_allowed_height = minimum_allowed_height self.forcing_terms.append(gravity) self.forcing_terms.append(manning_friction) 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 "level"' assert self.conserved_quantities[0] == 'level', msg msg = 'Second conserved quantity must be "xmomentum"' assert self.conserved_quantities[1] == 'xmomentum', msg msg = 'Third conserved quantity must be "ymomentum"' assert self.conserved_quantities[2] == 'ymomentum', msg #Check that levels are >= bed elevation from Numeric import alltrue, greater_equal level = self.quantities['level'] bed = self.quantities['elevation'] msg = 'All water levels must be greater than the bed elevation' assert alltrue( greater_equal( level.vertex_values, bed.vertex_values )), msg assert alltrue( greater_equal( level.edge_values, bed.edge_values )), msg assert alltrue( greater_equal( level.centroid_values, bed.centroid_values )), msg def compute_fluxes(self): #Call correct module function #(either from this module or C-extension) compute_fluxes(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) #################################### # Flux computation def flux_function(normal, ql, qr, zl, zr): """Compute fluxes between volumes for the shallow water wave equation cast in terms of w = h+z using the 'central scheme' as described in Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'. Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740. The implemented formula is given in equation (3.15) on page 714 Conserved quantities w, uh, vh are stored as elements 0, 1 and 2 in the numerical vectors ql an qr. Bed elevations zl and zr. #FIXME: Remove those and pass in height directly #(unless we'll include a bed elevation discontinuity later) """ from config import g, epsilon from math import sqrt from Numeric import array from util import rotate #Align momentums with x-axis q_left = rotate(ql, normal, direction = 1) q_right = rotate(qr, normal, direction = 1) z = (zl+zr)/2 #Take average of field values #FIXME: Maybe allow discontinuity later w_left = q_left[0] #w=h+z h_left = w_left-z uh_left = q_left[1] if h_left < epsilon: u_left = 0.0 #Could have been negative h_left = 0.0 else: u_left = uh_left/h_left w_right = q_right[0] #w=h+z h_right = w_right-z uh_right = q_right[1] if h_right < epsilon: u_right = 0.0 #Could have been negative h_right = 0.0 else: u_right = uh_right/h_right vh_left = q_left[2] vh_right = q_right[2] soundspeed_left = sqrt(g*h_left) soundspeed_right = sqrt(g*h_right) #Maximal wave speed s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0) #Minimal wave speed s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0) #Flux computation flux_left = array([u_left*h_left, u_left*uh_left + 0.5*g*h_left**2, u_left*vh_left]) flux_right = array([u_right*h_right, u_right*uh_right + 0.5*g*h_right**2, u_right*vh_right]) denom = s_max-s_min if denom == 0.0: edgeflux = array([0.0, 0.0, 0.0]) max_speed = 0.0 else: edgeflux = (s_max*flux_left - s_min*flux_right)/denom edgeflux += s_max*s_min*(q_right-q_left)/denom edgeflux = rotate(edgeflux, normal, direction=-1) max_speed = max(abs(s_max), abs(s_min)) return edgeflux, max_speed def compute_fluxes(domain): """Compute all fluxes and the timestep suitable for all volumes in domain. Compute total flux for each conserved quantity using "flux_function" Fluxes across each edge are scaled by edgelengths and summed up Resulting flux is then scaled by area and stored in domain.explicit_update The maximal allowable speed computed by the flux_function for each volume is converted to a timestep that must not be exceeded. The minimum of those is computed as the next overall timestep. Post conditions: domain.explicit_update is reset to computed flux values domain.timestep is set to the largest step satisfying all volumes. """ import sys from Numeric import zeros, Float from config import max_timestep N = domain.number_of_elements neighbours = domain.neighbours neighbour_edges = domain.neighbour_edges normals = domain.normals areas = domain.areas radii = domain.radii edgelengths = domain.edgelengths timestep = max_timestep #FIXME: Get rid of this #Shortcuts Level = domain.quantities['level'] Xmom = domain.quantities['xmomentum'] Ymom = domain.quantities['ymomentum'] Bed = domain.quantities['elevation'] #Arrays level = Level.edge_values xmom = Xmom.edge_values ymom = Ymom.edge_values bed = Bed.edge_values level_bdry = Level.boundary_values xmom_bdry = Xmom.boundary_values ymom_bdry = Ymom.boundary_values flux = zeros(3, Float) #Work array for summing up fluxes #Loop for k in range(N): optimal_timestep = float(sys.maxint) flux[:] = 0. #Reset work array for i in range(3): #Quantities inside volume facing neighbour i ql = [level[k, i], xmom[k, i], ymom[k, i]] zl = bed[k, i] #Quantities at neighbour on nearest face n = neighbours[k,i] if n < 0: m = -n-1 #Convert neg flag to index qr = [level_bdry[m], xmom_bdry[m], ymom_bdry[m]] zr = zl #Extend bed elevation to boundary else: m = neighbour_edges[k,i] qr = [level[n, m], xmom[n, m], ymom[n, m]] zr = bed[n, m] #Outward pointing normal vector normal = normals[k, 2*i:2*i+2] #Flux computation using provided function edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr) flux -= edgeflux * edgelengths[k,i] #Update optimal_timestep try: optimal_timestep = min(optimal_timestep, radii[k]/max_speed) except ZeroDivisionError: pass #Normalise by area and store for when all conserved #quantities get updated flux /= areas[k] Level.explicit_update[k] = flux[0] Xmom.explicit_update[k] = flux[1] Ymom.explicit_update[k] = flux[2] timestep = min(timestep, optimal_timestep) domain.timestep = timestep #################################### # Module functions for gradient limiting (distribute_to_vertices_and_edges) def distribute_to_vertices_and_edges(domain): ##print 'Calling distrib within SW' if domain.order == 1: #FIXME: This should be cleaned up, but we try to follow #pyvolution 2 strictly for now protect_against_negative_heights_centroid(domain) protect_against_infinitesimal_heights_centroid(domain) extrapolate_first_order(domain) elif domain.order == 2: #protect_against_negative_heights_centroid(domain) extrapolate_second_order(domain) else: raise 'Unknown order' #Compute edge values for name in domain.conserved_quantities: Q = domain.quantities[name] Q.interpolate_from_vertices_to_edges() def protect_against_infinitesimal_heights_centroid(domain): """Adjust height and momentum at centroid if height is less than minimum allowed height """ #FIXME: Used only in first order #Water levels at centroids wc = domain.quantities['level'].centroid_values #Bed elevations at centroids zc = domain.quantities['elevation'].centroid_values #Water depths at centroids hc = wc - zc #Momentums at centroids xmomc = domain.quantities['xmomentum'].centroid_values ymomc = domain.quantities['ymomentum'].centroid_values for k in range(domain.number_of_elements): #Protect against infinitesimal heights and high velocities if hc[k] < domain.minimum_allowed_height: #Control level and height if hc[k] < 0.0: wc[k] = zc[k]; hc[k] = 0.0 #Control momentum xmomc[k] = ymomc[k] = 0.0 def protect_against_negative_heights_centroid(domain): """Adjust height and momentum at centroid if height is less than zero """ #Water levels at centroids wc = domain.quantities['level'].centroid_values #Bed elevations at centroids zc = domain.quantities['elevation'].centroid_values #Water depths at centroids hc = wc - zc #Momentums at centroids xmomc = domain.quantities['xmomentum'].centroid_values ymomc = domain.quantities['ymomentum'].centroid_values for k in range(domain.number_of_elements): #Protect against infinitesimal heights and high velocities if hc[k] < 0.0: #Control level and height wc[k] = zc[k]; hc[k] = 0.0 #Control momentum xmomc[k] = ymomc[k] = 0.0 def extrapolate_first_order(domain): """First order extrapolator function, 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 (see comment in code) In addition, momemtums get distributed as constant values. Precondition: All quantities defined at centroids and bed elevation defined at vertices. Postcondition Conserved quantities defined at vertices """ #Update conserved quantities using straight first order for name in domain.conserved_quantities: Q = domain.quantities[name] Q.extrapolate_first_order() #Water levels at centroids wc = domain.quantities['level'].centroid_values #Bed elevations at centroids zc = domain.quantities['elevation'].centroid_values #Water depths at centroids hc = wc - zc #Water levels at vertices wv = domain.quantities['level'].vertex_values #Bed elevations at vertices zv = domain.quantities['elevation'].vertex_values #Water depths at vertices hv = wv-zv #Computed weighted balance between constant levels and and #levels parallel to the bed elevation. for k in range(domain.number_of_elements): #Compute maximal variation in bed elevation z_range = max(abs(zv[k,0]-zc[k]), abs(zv[k,1]-zc[k]), abs(zv[k,2]-zc[k])) #Weighted balance between stage parallel to bed elevation #(wvi = zvi + hc) and constant stage (wvi = wc = zc+hc) #where i=0,1,2 denotes the vertex ids # #It follows that # wvi = (1-alpha)*(zvi+hc) + alpha*(zc+hc) = # (1-alpha)*zvi + alpha*zc + hc = # zvi + hc + alpha*(zc-zvi) # #where alpha in [0,1] and defined as the ratio between hc and #the maximal difference from zc to zv0, zv1 and zv2 # #Mathematically the following can be continued on using hc as # wvi = # zvi + hc + alpha*(zc+hc-zvi-hc) = # zvi + hc + alpha*(hvi-hc) #since hvi = zc+hc-zvi in the constant case if z_range > 0.0: alpha = min(hc[k]/z_range, 1.0) else: alpha = 1.0 #Update water levels for i in range(3): #FIXME: Use the original first-order one first, then switch wv[k,i] = zv[k,i] + hc[k] + alpha*(zc[k]-zv[k,i]) #wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k]) #FIXME: What about alpha weighting of momentum?? def extrapolate_second_order(domain): """Second order limiter function, 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 (see comment in code) A weighted average between shallow and deep cases is as in the first order case. In addition, all conserved quantities get distributed as per a piecewise linear function. 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 """ #FIXME: first and second order might merge #Update conserved quantities using straight second order for name in domain.conserved_quantities: Q = domain.quantities[name] Q.extrapolate_second_order() #print 'y1', Q.vertex_values[1,:] #OK #FIXME - like pyvolution 2 ....................... protect_against_negative_heights_centroid(domain) #print 'y1', Q.vertex_values[1,:] #OK for name in domain.conserved_quantities: Q = domain.quantities[name] Q.limit() #print 'y1', Q.vertex_values[1,:] #OK #Water levels at centroids wc = domain.quantities['level'].centroid_values #Bed elevations at centroids zc = domain.quantities['elevation'].centroid_values #Water depths at centroids hc = wc - zc #Water levels at vertices wv = domain.quantities['level'].vertex_values #Bed elevations at vertices zv = domain.quantities['elevation'].vertex_values #Water depths at vertices hv = wv-zv #Computed linear combination between constant levels and and #levels parallel to the bed elevation. 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 # z_range = 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 shallow #first order scheme and alpha==1 means using the limited #second order scheme for the stage w #If hmin < 0 then alpha=0 reverting to first order. if z_range > 0.0: alpha = max( min( 2*hmin/z_range, 1.0), 0.0) else: alpha = 1.0 ##if k==1: print 'alpha', alpha #OK #Update water levels # (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) = # zvi + hc + alpha*(hvi - hc) if alpha < 1: for i in range(3): wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k]) #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-alpha)*xmomc[k] + alpha*xmomv[k,:]; ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]; #print 'y1', Q.vertex_values[1,:] #OK #Finally, protect against infinitesimal heights and high speeds #Water depths at vertices hv = wv-zv hc = wc-zc for k in range(domain.number_of_elements): hmax = max(hv[k,:]) if hmax < domain.minimum_allowed_height: #Reset negative heights to bed elevation if hc[k] < 0.0: wc[k] = zc[k] for i in range(3): if hv[k,i] < 0.0: wv[k,i] = zv[k,i] ############################################### #Boundary - 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 self.domain = domain 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.domain.get_conserved_quantities(vol_id, edge = edge_id) normal = self.domain.get_normal(vol_id, edge_id) r = rotate(q, normal, direction = 1) r[1] = -r[1] q = rotate(r, normal, direction = -1) return q ######################### #Standard forcing terms: # def gravity(domain): """Implement forcing function for bed slope working with consecutive data structures of class Volume """ from config import g from util import gradient from Numeric import zeros, Float, array, sum Level = domain.quantities['level'] Xmom = domain.quantities['xmomentum'] Ymom = domain.quantities['ymomentum'] Elevation = domain.quantities['elevation'] h = Level.edge_values - Elevation.edge_values V = domain.get_vertex_coordinates() for k in range(domain.number_of_elements): avg_h = sum( h[k,:] )/3 #Compute bed slope x0, y0, x1, y1, x2, y2 = V[k,:] z0, z1, z2 = Elevation.vertex_values[k,:] zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2) #Update momentum Xmom.explicit_update[k] += -g*zx*avg_h Ymom.explicit_update[k] += -g*zy*avg_h def manning_friction(domain): """Apply (Manning) friction to water momentum """ import Numeric from math import sqrt from config import g, minimum_allowed_height Level = domain.quantities['level'] Xmom = domain.quantities['xmomentum'] Ymom = domain.quantities['ymomentum'] Friction = domain.quantities['friction'] for k in range(domain.number_of_elements): w = Level.centroid_values[k] uh = Xmom.centroid_values[k] vh = Ymom.centroid_values[k] manning = Friction.centroid_values[k] if w >= minimum_allowed_height: S = -g * manning**2 * sqrt((uh**2 + vh**2)) S /= w**(7.0/3) #Update momentum Xmom.explicit_update[k] += S*uh Ymom.explicit_update[k] += S*vh ########################### ########################### #Geometries #FIXME: Rethink this way of creating values. class Weir: """Set a bathymetry for weir with a hole and a downstream gutter x,y are assumed to be in the unit square """ def __init__(self, stage): self.inflow_stage = stage def __call__(self, x, y): from Numeric import zeros, Float from math import sqrt N = len(x) assert N == len(y) z = zeros(N, Float) for i in range(N): z[i] = -x[i]/2 #General slope #Flattish bit to the left if x[i] < 0.3: z[i] = -x[i]/10 #Weir if x[i] >= 0.3 and x[i] < 0.4: z[i] = -x[i]+0.9 #Dip x0 = 0.6 #depth = -1.3 depth = -1.0 #plateaux = -0.9 plateaux = -0.6 if y[i] < 0.7: if x[i] > x0 and x[i] < 0.9: z[i] = depth #RHS plateaux if x[i] >= 0.9: z[i] = plateaux elif y[i] >= 0.7 and y[i] < 1.5: #Restrict and deepen if x[i] >= x0 and x[i] < 0.8: z[i] = depth-(y[i]/3-0.3) #z[i] = depth-y[i]/5 #z[i] = depth elif x[i] >= 0.8: #RHS plateaux z[i] = plateaux elif y[i] >= 1.5: if x[i] >= x0 and x[i] < 0.8 + (y[i]-1.5)/1.2: #Widen up and stay at constant depth z[i] = depth-1.5/5 elif x[i] >= 0.8 + (y[i]-1.5)/1.2: #RHS plateaux z[i] = plateaux #Hole in weir (slightly higher than inflow condition) if x[i] >= 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4: z[i] = -x[i]+self.inflow_stage + 0.02 #Channel behind weir x0 = 0.5 if x[i] >= 0.4 and x[i] < x0 and y[i] > 0.2 and y[i] < 0.4: z[i] = -x[i]+self.inflow_stage + 0.02 if x[i] >= x0 and x[i] < 0.6 and y[i] > 0.2 and y[i] < 0.4: #Flatten it out towards the end z[i] = -x0+self.inflow_stage + 0.02 + (x0-x[i])/5 #Hole to the east x0 = 1.1; y0 = 0.35 #if x[i] < -0.2 and y < 0.5: if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2: z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-1.0 #Tiny channel draining hole if x[i] >= 1.14 and x[i] < 1.2 and y[i] >= 0.4 and y[i] < 0.6: z[i] = -0.9 #North south if x[i] >= 0.9 and x[i] < 1.18 and y[i] >= 0.58 and y[i] < 0.65: z[i] = -1.0 + (x[i]-0.9)/3 #East west #Stuff not in use #Upward slope at inlet to the north west #if x[i] < 0.0: # and y[i] > 0.5: # #z[i] = -y[i]+0.5 #-x[i]/2 # z[i] = x[i]/4 - y[i]**2 + 0.5 #Hole to the west #x0 = -0.4; y0 = 0.35 # center #if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2: # z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-0.2 return z/2 class Weir_simple: """Set a bathymetry for weir with a hole and a downstream gutter x,y are assumed to be in the unit square """ def __init__(self, stage): self.inflow_stage = stage def __call__(self, x, y): from Numeric import zeros, Float N = len(x) assert N == len(y) z = zeros(N, Float) for i in range(N): z[i] = -x[i] #General slope #Flat bit to the left if x[i] < 0.3: z[i] = -x[i]/10 #General slope #Weir if x[i] > 0.3 and x[i] < 0.4: z[i] = -x[i]+0.9 #Dip if x[i] > 0.6 and x[i] < 0.9: z[i] = -x[i]-0.5 #-y[i]/5 #Hole in weir (slightly higher than inflow condition) if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4: z[i] = -x[i]+self.inflow_stage + 0.05 return z/2 class Constant_height: """Set an initial condition with constant water height, e.g stage s = z+h """ def __init__(self, W, h): self.W = W self.h = h def __call__(self, x, y): if self.W is None: from Numeric import ones, Float return self.h*ones(len(x), Float) else: return self.W(x,y) + self.h ############################################## #Initialise module