# Changeset 4769

Ignore:
Timestamp:
Oct 30, 2007, 5:13:17 PM (16 years ago)
Message:

Retired all Python code which is superseded by C implementations.
This has made the code base much smaller and there is less confusion when people try to understand the algorithms.

Location:
anuga_core/source/anuga
Files:
7 edited

Unmodified
Removed
• ## anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

 r4768 #General input checks # General input checks L = [numeric, quantity, function, geospatial_data, points, filename] msg = 'Exactly one of the arguments '+\ #Determine which 'set_values_from_...' to use # Determine which 'set_values_from_...' to use if numeric is not None: def update(quantity, timestep): """Update centroid values based on values stored in explicit_update and semi_implicit_update as well as given timestep Function implementing forcing terms must take on argument which is the domain and they must update either explicit or implicit updates, e,g,: def gravity(domain): .... domain.quantities['xmomentum'].explicit_update = ... domain.quantities['ymomentum'].explicit_update = ... Explicit terms must have the form G(q, t) and explicit scheme is q^{(n+1}) = q^{(n)} + delta_t G(q^{n}, n delta_t) Semi implicit forcing terms are assumed to have the form G(q, t) = H(q, t) q and the semi implicit scheme will then be q^{(n+1}) = q^{(n)} + delta_t H(q^{n}, n delta_t) q^{(n+1}) """ from Numeric import sum, equal, ones, exp, Float N = quantity.centroid_values.shape[0] # Divide H by conserved quantity to obtain G (see docstring above) for k in range(N): x = quantity.centroid_values[k] if x == 0.0: #FIXME: Is this right quantity.semi_implicit_update[k] = 0.0 else: quantity.semi_implicit_update[k] /= x # Semi implicit updates denominator = ones(N, Float)-timestep*quantity.semi_implicit_update if sum(less(denominator, 1.0)) > 0.0: msg = 'denominator < 1.0 in semi implicit update. Call Stephen :-)' raise msg if sum(equal(denominator, 0.0)) > 0.0: msg = 'Zero division in semi implicit update. Call Stephen :-)' raise msg else: #Update conserved_quantities from semi implicit updates quantity.centroid_values /= denominator #    quantity.centroid_values = exp(timestep*quantity.semi_implicit_update)*quantity.centroid_values #Explicit updates quantity.centroid_values += timestep*quantity.explicit_update def interpolate_from_vertices_to_edges(quantity): """Compute edge values from vertex values using linear interpolation """ for k in range(quantity.vertex_values.shape[0]): q0 = quantity.vertex_values[k, 0] q1 = quantity.vertex_values[k, 1] q2 = quantity.vertex_values[k, 2] quantity.edge_values[k, 0] = 0.5*(q1+q2) quantity.edge_values[k, 1] = 0.5*(q0+q2) quantity.edge_values[k, 2] = 0.5*(q0+q1) def backup_centroid_values(quantity): """Copy centroid values to backup array""" qc = quantity.centroid_values qb = quantity.centroid_backup_values #Check each triangle for k in range(len(quantity.domain)): qb[k] = qc[k] def saxpy_centroid_values(quantity,a,b): """saxpy operation between centroid value and backup""" qc = quantity.centroid_values qb = quantity.centroid_backup_values #Check each triangle for k in range(len(quantity.domain)): qc[k] = a*qc[k]+b*qb[k] def compute_gradients(quantity): """Compute gradients of triangle surfaces defined by centroids of neighbouring volumes. If one edge is on the boundary, use own centroid as neighbour centroid. If two or more are on the boundary, fall back to first order scheme. """ from Numeric import zeros, Float from utilitites.numerical_tools import gradient centroid_coordinates = quantity.domain.centroid_coordinates surrogate_neighbours = quantity.domain.surrogate_neighbours centroid_values = quantity.centroid_values number_of_boundaries = quantity.domain.number_of_boundaries N = centroid_values.shape[0] a = zeros(N, Float) b = zeros(N, Float) for k in range(N): if number_of_boundaries[k] < 2: #Two or three true neighbours #Get indices of neighbours (or self when used as surrogate) k0, k1, k2 = surrogate_neighbours[k,:] #Get data q0 = centroid_values[k0] q1 = centroid_values[k1] q2 = centroid_values[k2] x0, y0 = centroid_coordinates[k0] #V0 centroid x1, y1 = centroid_coordinates[k1] #V1 centroid x2, y2 = centroid_coordinates[k2] #V2 centroid #Gradient a[k], b[k] = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2) elif number_of_boundaries[k] == 2: #One true neighbour #Get index of the one neighbour for k0 in surrogate_neighbours[k,:]: if k0 != k: break assert k0 != k k1 = k  #self #Get data q0 = centroid_values[k0] q1 = centroid_values[k1] x0, y0 = centroid_coordinates[k0] #V0 centroid x1, y1 = centroid_coordinates[k1] #V1 centroid #Gradient a[k], b[k] = gradient2(x0, y0, x1, y1, q0, q1) else: #No true neighbours - #Fall back to first order scheme pass return a, b def limit(quantity): """Limit slopes for each volume to eliminate artificial variance introduced by e.g. second order extrapolator This is an unsophisticated limiter as it does not take into account dependencies among quantities. precondition: vertex values are estimated from gradient postcondition: vertex values are updated """ from Numeric import zeros, Float N = quantity.domain.number_of_nodes beta_w = quantity.domain.beta_w qc = quantity.centroid_values qv = quantity.vertex_values #Find min and max of this and neighbour's centroid values qmax = zeros(qc.shape, Float) qmin = zeros(qc.shape, Float) for k in range(N): qmax[k] = qc[k] qmin[k] = qc[k] for i in range(3): n = quantity.domain.neighbours[k,i] if n >= 0: qn = qc[n] #Neighbour's centroid value qmin[k] = min(qmin[k], qn) qmax[k] = max(qmax[k], qn) qmax[k] = min(qmax[k], 2.0*qc[k]) qmin[k] = max(qmin[k], 0.5*qc[k]) #Diffences between centroids and maxima/minima dqmax = qmax - qc dqmin = qmin - qc #Deltas between vertex and centroid values dq = zeros(qv.shape, Float) for i in range(3): dq[:,i] = qv[:,i] - qc #Phi limiter for k in range(N): #Find the gradient limiter (phi) across vertices phi = 1.0 for i in range(3): r = 1.0 if (dq[k,i] > 0): r = dqmax[k]/dq[k,i] if (dq[k,i] < 0): r = dqmin[k]/dq[k,i] phi = min( min(r*beta_w, 1), phi ) #Then update using phi limiter for i in range(3): qv[k,i] = qc[k] + phi*dq[k,i] from anuga.utilities import compile if compile.can_use_C_extension('quantity_ext.c'): #Replace python version with c implementations if compile.can_use_C_extension('quantity_ext.c'): # Underlying C implementations can be accessed from quantity_ext import average_vertex_values, backup_centroid_values, \ from quantity_ext import compute_gradients, limit,\ extrapolate_second_order, interpolate_from_vertices_to_edges, update extrapolate_second_order, interpolate_from_vertices_to_edges, update else: msg = 'C implementations could not be accessed by %s.\n ' %__file__ msg += 'Make sure compile_all.py has been run as described in ' msg += 'the ANUGA installation guide.' raise Exception, msg
• ## anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

 r4728 // Gateways to Python PyObject *update(PyObject *self, PyObject *args) { // FIXME (Ole): It would be great to turn this text into a Python DOC string /*"""Update centroid values based on values stored in explicit_update and semi_implicit_update as well as given timestep Function implementing forcing terms must take on argument which is the domain and they must update either explicit or implicit updates, e,g,: def gravity(domain): .... domain.quantities['xmomentum'].explicit_update = ... domain.quantities['ymomentum'].explicit_update = ... Explicit terms must have the form G(q, t) and explicit scheme is q^{(n+1}) = q^{(n)} + delta_t G(q^{n}, n delta_t) Semi implicit forcing terms are assumed to have the form G(q, t) = H(q, t) q and the semi implicit scheme will then be q^{(n+1}) = q^{(n)} + delta_t H(q^{n}, n delta_t) q^{(n+1}) */ PyObject *quantity; PyObject *interpolate_from_vertices_to_edges(PyObject *self, PyObject *args) { // //Compute edge values from vertex values using linear interpolation // PyObject *quantity; PyArrayObject *vertex_values, *edge_values; PyObject *compute_gradients(PyObject *self, PyObject *args) { //"""Compute gradients of triangle surfaces defined by centroids of //neighbouring volumes. //If one edge is on the boundary, use own centroid as neighbour centroid. //If two or more are on the boundary, fall back to first order scheme. //""" PyObject *quantity, *domain, *R; PyObject *limit(PyObject *self, PyObject *args) { //Limit slopes for each volume to eliminate artificial variance //introduced by e.g. second order extrapolator //This is an unsophisticated limiter as it does not take into //account dependencies among quantities. //precondition: //  vertex values are estimated from gradient //postcondition: //  vertex values are updated // PyObject *quantity, *domain, *Tmp;
• ## anuga_core/source/anuga/shallow_water/shallow_water_domain.py

 r4754 """ #Subversion keywords: # Subversion keywords: # #\$LastChangedDate\$ #\$LastChangedRevision\$ #\$LastChangedBy\$ # \$LastChangedDate\$ # \$LastChangedRevision\$ # \$LastChangedBy\$ from Numeric import zeros, ones, Float, array, sum, size from anuga.config import optimise_dry_cells #Shallow water domain #--------------------- # Shallow water domain #--------------------- class Domain(Generic_Domain): number_of_full_triangles=number_of_full_triangles) #self.minimum_allowed_height = minimum_allowed_height #self.H0 = minimum_allowed_height self.set_minimum_allowed_height(minimum_allowed_height) self.optimise_dry_cells = optimise_dry_cells self.flux_function = flux_function_central #self.flux_function = flux_function_kinetic self.forcing_terms.append(manning_friction) self.forcing_terms.append(manning_friction_implicit) self.forcing_terms.append(gravity) def distribute_to_vertices_and_edges(self): #Call correct module function #(either from this module or C-extension) # Call correct module function # (either from this module or C-extension) distribute_to_vertices_and_edges(self) #FIXME: Under construction #     def set_defaults(self): #         """Set default values for uninitialised quantities. #         This is specific to the shallow water wave equation #         Defaults for 'elevation', 'friction', 'xmomentum' and 'ymomentum' #         are 0.0. Default for 'stage' is whatever the value of 'elevation'. #         """ #         for name in self.other_quantities + self.conserved_quantities: #             print name #             print self.quantities.keys() #             if not self.quantities.has_key(name): #                 if name == 'stage': #                     if self.quantities.has_key('elevation'): #                         z = self.quantities['elevation'].vertex_values #                         self.set_quantity(name, z) #                     else: #                         self.set_quantity(name, 0.0) #                 else: #                     self.set_quantity(name, 0.0) #         #Lift negative heights up #         #z = self.quantities['elevation'].vertex_values #         #w = self.quantities['stage'].vertex_values #         #h = w-z #         #for k in range(h.shape[0]): #         #    for i in range(3): #         #        if h[k, i] < 0.0: #         #            w[k, i] = z[k, i] #         #self.quantities['stage'].interpolate() """ #Call check integrity here rather than from user scripts #self.check_integrity() # Call check integrity here rather than from user scripts # self.check_integrity() msg = 'Parameter beta_h must be in the interval [0, 1[' #Initial update of vertex and edge values before any STORAGE #and or visualisation #This is done again in the initialisation of the Generic_Domain #evolve loop but we do it here to ensure the values are ok for storage # Initial update of vertex and edge values before any STORAGE # and or visualisation # This is done again in the initialisation of the Generic_Domain # evolve loop but we do it here to ensure the values are ok for storage self.distribute_to_vertices_and_edges() if self.store is True and self.time == 0.0: self.initialise_storage() #print 'Storing results in ' + self.writer.filename # print 'Storing results in ' + self.writer.filename else: pass #print 'Results will not be stored.' #print 'To store results set domain.store = True' #FIXME: Diagnostic output should be controlled by # print 'Results will not be stored.' # print 'To store results set domain.store = True' # FIXME: Diagnostic output should be controlled by # a 'verbose' flag living in domain (or in a parent class) #Call basic machinery from parent class # Call basic machinery from parent class for t in Generic_Domain.evolve(self, yieldstep=yieldstep, skip_initial_step=skip_initial_step): #Store model data, e.g. for subsequent visualisation # Store model data, e.g. for subsequent visualisation if self.store is True: self.store_timestep(self.quantities_to_be_stored) #FIXME: Could maybe be taken from specified list #of 'store every step' quantities #Pass control on to outer loop for more specific actions # FIXME: Could maybe be taken from specified list # of 'store every step' quantities # Pass control on to outer loop for more specific actions yield(t) def initialise_storage(self): from anuga.shallow_water.data_manager import get_dataobject #Initialise writer # Initialise writer self.writer = get_dataobject(self, mode = 'w') #Store vertices and connectivity # Store vertices and connectivity self.writer.store_connectivity() #=============== 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 """ 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 #################################### #=============== End of class Shallow Water Domain =============================== #----------------- # Flux computation # FIXME (Ole): Delete Python versions of code that is # always used in C-extensions - they get out of sync anyway def flux_function_central(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 and qr. Bed elevations zl and zr. """ from math import sqrt #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 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 #FIXME(Ole): Why is it again that we don't #use uh_left and uh_right directly in the first entries? 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 erfcc(x): from math import fabs, exp z=fabs(x) t=1.0/(1.0+0.5*z) result=t*exp(-z*z-1.26551223+t*(1.00002368+t*(.37409196+ t*(.09678418+t*(-.18628806+t*(.27886807+t*(-1.13520398+ t*(1.48851587+t*(-.82215223+t*.17087277))))))))) if x < 0.0: result = 2.0-result return result def flux_function_kinetic(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 Zhang et. al., Advances in Water Resources, 26(6), 2003, 635-647. 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. """ from math import sqrt from Numeric import array #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 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 F_left  = 0.0 F_right = 0.0 from math import sqrt, pi, exp if h_left > 0.0: F_left = u_left/sqrt(g*h_left) if h_right > 0.0: F_right = u_right/sqrt(g*h_right) edgeflux = array([0.0, 0.0, 0.0]) edgeflux[0] = h_left*u_left/2.0*erfcc(-F_left) +  \ h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left**2)) + \ h_right*u_right/2.0*erfcc(F_right) -  \ h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right**2)) edgeflux[1] = (h_left*u_left**2 + g/2.0*h_left**2)/2.0*erfcc(-F_left) + \ u_left*h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left**2)) + \ (h_right*u_right**2 + g/2.0*h_right**2)/2.0*erfcc(F_right) -  \ u_right*h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right**2)) edgeflux[2] = vh_left*u_left/2.0*erfcc(-F_left) + \ vh_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left**2)) + \ vh_right*u_right/2.0*erfcc(F_right) - \ vh_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right**2)) edgeflux = rotate(edgeflux, normal, direction=-1) max_speed = max(abs(s_max), abs(s_min)) return edgeflux, max_speed #----------------- def compute_fluxes(domain): domain.explicit_update is reset to computed flux values domain.timestep is set to the largest step satisfying all volumes. This wrapper calls the underlying C version of compute fluxes """ N = len(domain) # number_of_triangles #Shortcuts # Shortcuts Stage = domain.quantities['stage'] Xmom = domain.quantities['xmomentum'] Bed = domain.quantities['elevation'] #Arrays stage = Stage.edge_values xmom =  Xmom.edge_values ymom =  Ymom.edge_values bed =   Bed.edge_values stage_bdry = Stage.boundary_values xmom_bdry =  Xmom.boundary_values ymom_bdry =  Ymom.boundary_values flux = zeros(3, Float) #Work array for summing up fluxes #Loop timestep = float(sys.maxint) for k in range(N): flux[:] = 0.  #Reset work array for i in range(3): #Quantities inside volume facing neighbour i ql = [stage[k, i], xmom[k, i], ymom[k, i]] zl = bed[k, i] #Quantities at neighbour on nearest face n = domain.neighbours[k,i] if n < 0: m = -n-1 #Convert negative flag to index qr = [stage_bdry[m], xmom_bdry[m], ymom_bdry[m]] zr = zl #Extend bed elevation to boundary else: m = domain.neighbour_edges[k,i] qr = [stage[n, m], xmom[n, m], ymom[n, m]] zr = bed[n, m] #Outward pointing normal vector normal = domain.normals[k, 2*i:2*i+2] #Flux computation using provided function edgeflux, max_speed = domain.flux_function(normal, ql, qr, zl, zr) flux -= edgeflux * domain.edgelengths[k,i] #Update optimal_timestep on full cells if  domain.tri_full_flag[k] == 1: try: timestep = min(timestep, 0.5*domain.radii[k]/max_speed) except ZeroDivisionError: pass #Normalise by area and store for when all conserved #quantities get updated flux /= domain.areas[k] Stage.explicit_update[k] = flux[0] Xmom.explicit_update[k] = flux[1] Ymom.explicit_update[k] = flux[2] domain.max_speed[k] = max_speed domain.flux_timestep = timestep #MH090605 The following method belongs to the shallow_water domain class #see comments in the corresponding method in shallow_water_ext.c def extrapolate_second_order_sw_c(domain): from shallow_water_ext import\ compute_fluxes_ext_central as compute_fluxes_ext flux_timestep = compute_fluxes_ext(timestep, domain.epsilon, domain.H0, domain.g, domain.neighbours, domain.neighbour_edges, domain.normals, domain.edgelengths, domain.radii, domain.areas, domain.tri_full_flag, Stage.edge_values, Xmom.edge_values, Ymom.edge_values, Bed.edge_values, Stage.boundary_values, Xmom.boundary_values, Ymom.boundary_values, Stage.explicit_update, Xmom.explicit_update, Ymom.explicit_update, domain.already_computed_flux, domain.max_speed, int(domain.optimise_dry_cells)) domain.flux_timestep = flux_timestep #--------------------------------------- # Module functions for gradient limiting #--------------------------------------- # MH090605 The following method belongs to the shallow_water domain class # see comments in the corresponding method in shallow_water_ext.c def extrapolate_second_order_sw(domain): """Wrapper calling C version of extrapolate_second_order_sw """ Elevation = domain.quantities['elevation'] from shallow_water_ext import extrapolate_second_order_sw extrapolate_second_order_sw(domain, domain.surrogate_neighbours, domain.number_of_boundaries, domain.centroid_coordinates, Stage.centroid_values, Xmom.centroid_values, Ymom.centroid_values, Elevation.centroid_values, domain.vertex_coordinates, Stage.vertex_values, Xmom.vertex_values, Ymom.vertex_values, Elevation.vertex_values, int(domain.optimise_dry_cells)) def compute_fluxes_c(domain): """Wrapper calling C version of compute fluxes """ import sys N = len(domain) # number_of_triangles #Shortcuts Stage = domain.quantities['stage'] Xmom = domain.quantities['xmomentum'] Ymom = domain.quantities['ymomentum'] Bed = domain.quantities['elevation'] timestep = float(sys.maxint) from shallow_water_ext import\ compute_fluxes_ext_central as compute_fluxes_ext domain.flux_timestep = compute_fluxes_ext(timestep, domain.epsilon, domain.H0, domain.g, domain.neighbours, domain.neighbour_edges, domain.normals, domain.edgelengths, domain.radii, domain.areas, domain.tri_full_flag, Stage.edge_values, Xmom.edge_values, Ymom.edge_values, Bed.edge_values, Stage.boundary_values, Xmom.boundary_values, Ymom.boundary_values, Stage.explicit_update, Xmom.explicit_update, Ymom.explicit_update, domain.already_computed_flux, domain.max_speed, int(domain.optimise_dry_cells)) #################################### # Module functions for gradient limiting (distribute_to_vertices_and_edges) from shallow_water_ext import extrapolate_second_order_sw as extrapol2 extrapol2(domain, domain.surrogate_neighbours, domain.number_of_boundaries, domain.centroid_coordinates, Stage.centroid_values, Xmom.centroid_values, Ymom.centroid_values, Elevation.centroid_values, domain.vertex_coordinates, Stage.vertex_values, Xmom.vertex_values, Ymom.vertex_values, Elevation.vertex_values, int(domain.optimise_dry_cells)) def distribute_to_vertices_and_edges(domain): from anuga.config import optimised_gradient_limiter #Remove very thin layers of water # Remove very thin layers of water protect_against_infinitesimal_and_negative_heights(domain) #Extrapolate all conserved quantities # Extrapolate all conserved quantities if optimised_gradient_limiter: #MH090605 if second order, #perform the extrapolation and limiting on #all of the conserved quantities # MH090605 if second order, # perform the extrapolation and limiting on # all of the conserved quantities if (domain._order_ == 1): raise 'Unknown order' else: #old code: # Old code: for name in domain.conserved_quantities: Q = domain.quantities[name] """ #Shortcuts # 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 #Update #FIXME: Modify according to c-version - or discard altogether. for k in range(len(domain)): if hc[k] < domain.minimum_allowed_height: #Control stage if hc[k] < domain.epsilon: wc[k] = zc[k] # Contain 'lost mass' error #Control momentum xmomc[k] = ymomc[k] = 0.0 def protect_against_infinitesimal_and_negative_heights_c(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 from shallow_water_ext import protect protect(domain.minimum_allowed_height, domain.maximum_allowed_speed, domain.epsilon, wc, zc, xmomc, ymomc) This limiter depends on two quantities (w,z) so it resides within this module rather than within quantity.py """ N = len(domain) Wrapper for c-extension """ N = len(domain) # number_of_triangles beta_h = domain.beta_h #Shortcuts # Shortcuts wc = domain.quantities['stage'].centroid_values zc = domain.quantities['elevation'].centroid_values 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): 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): 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): 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): hvbar[k,i] = hc[k] + phi*dh[k,i] return hvbar def h_limiter_c(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 Wrapper for c-extension """ N = len(domain) # number_of_triangles 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 #Call C-extension from shallow_water_ext import h_limiter_sw as h_limiter hvbar = h_limiter(domain, hc, hv) from shallow_water_ext import h_limiter_sw hvbar = h_limiter_sw(domain, hc, hv) return hvbar 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(len(domain)): #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 #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): wv[k,i] = zv[k,i] + (1-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-alpha)*xmomc[k] + alpha*xmomv[k,:] ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:] def balance_deep_and_shallow_c(domain): """Wrapper for C implementation Wrapper for C implementation """ # Omit updating xmomv # from shallow_water_ext import balance_deep_and_shallow from shallow_water_ext import balance_deep_and_shallow as balance_deep_and_shallow_c # 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 # Momentums at centroids ymomv = domain.quantities['ymomentum'].vertex_values #Limit h # Limit h if domain.beta_h > 0: hvbar = h_limiter(domain) balance_deep_and_shallow(domain, domain.beta_h, wc, zc, wv, zv, hvbar, xmomc, ymomc, xmomv, ymomv) balance_deep_and_shallow_c(domain, domain.beta_h, wc, zc, wv, zv, hvbar, xmomc, ymomc, xmomv, ymomv) else: # print 'Using first order h-limiter' # FIXME: Pass wc in for now - it will be ignored. #This is how one would make a first order h_limited value #as in the old balancer (pre 17 Feb 2005): # If we wish to hard wire this, one should modify the C-code #from Numeric import zeros, Float #hvbar = zeros( (len(wc), 3), Float) #for i in range(3): #    hvbar[:,i] = wc[:] - zc[:] balance_deep_and_shallow(domain, domain.beta_h, wc, zc, wv, zv, wc, xmomc, ymomc, xmomv, ymomv) ############################################### #Boundaries - specific to the shallow water wave equation # This is how one would make a first order h_limited value # as in the old balancer (pre 17 Feb 2005): #  If we wish to hard wire this, one should modify the C-code # from Numeric import zeros, Float # hvbar = zeros( (len(wc), 3), Float) # for i in range(3): #     hvbar[:,i] = wc[:] - zc[:] balance_deep_and_shallow_c(domain, domain.beta_h, wc, zc, wv, zv, wc, xmomc, ymomc, xmomv, ymomv) #------------------------------------------------------------------ # Boundary conditions - specific to the shallow water wave equation #------------------------------------------------------------------ class Reflective_boundary(Boundary): """Reflective boundary returns same conserved quantities as raise msg #Handy shorthands # Handy shorthands self.stage   = domain.quantities['stage'].edge_values self.xmom    = domain.quantities['xmomentum'].edge_values #FIXME: Consider this (taken from File_boundary) to allow #spatial variation #if vol_id is not None and edge_id is not None: #    i = self.boundary_indices[ vol_id, edge_id ] #    return self.F(t, point_id = i) #else: #    return self.F(t) # FIXME: Consider this (taken from File_boundary) to allow # spatial variation # if vol_id is not None and edge_id is not None: #     i = self.boundary_indices[ vol_id, edge_id ] #     return self.F(t, point_id = i) # else: #     return self.F(t) #FIXME: Consider this (taken from File_boundary) to allow #spatial variation #if vol_id is not None and edge_id is not None: #    i = self.boundary_indices[ vol_id, edge_id ] #    return self.F(t, point_id = i) #else: #    return self.F(t) # FIXME: Consider this (taken from File_boundary) to allow # spatial variation # if vol_id is not None and edge_id is not None: #     i = self.boundary_indices[ vol_id, edge_id ] #     return self.F(t, point_id = i) # else: #     return self.F(t) ######################### #Standard forcing terms: # #----------------------- # Standard forcing terms #----------------------- def gravity(domain): """Apply gravitational pull in the presence of bed slope Wrapper calls underlying C implementation """ ymom = domain.quantities['ymomentum'].explicit_update Stage = domain.quantities['stage'] Elevation = domain.quantities['elevation'] h = Stage.edge_values - Elevation.edge_values v = Elevation.vertex_values stage = domain.quantities['stage'] elevation = domain.quantities['elevation'] h = stage.centroid_values - elevation.centroid_values z = elevation.vertex_values x = domain.get_vertex_coordinates() g = domain.g for k in range(len(domain)): avg_h = sum( h[k,:] )/3 #Compute bed slope x0, y0, x1, y1, x2, y2 = x[k,:] z0, z1, z2 = v[k,:] zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2) #Update momentum xmom[k] += -g*zx*avg_h ymom[k] += -g*zy*avg_h def gravity_c(domain): """Wrapper calling C version """ xmom = domain.quantities['xmomentum'].explicit_update ymom = domain.quantities['ymomentum'].explicit_update stage = domain.quantities['stage'] elevation = domain.quantities['elevation'] h = stage.centroid_values - elevation.centroid_values z = elevation.vertex_values x = domain.get_vertex_coordinates() g = domain.g from shallow_water_ext import gravity gravity(g, h, z, x, xmom, ymom) #, 1.0e-6) def manning_friction(domain): """Apply (Manning) friction to water momentum (Python version) """ from math import sqrt from shallow_water_ext import gravity as gravity_c gravity_c(g, h, z, x, xmom, ymom) #, 1.0e-6) def manning_friction_implicit(domain): """Apply (Manning) friction to water momentum Wrapper for c version """ #print 'Implicit friction' xmom = domain.quantities['xmomentum'] ymom = domain.quantities['ymomentum'] 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 uh = xmom.centroid_values vh = ymom.centroid_values eta = domain.quantities['friction'].centroid_values xmom_update = domain.quantities['xmomentum'].semi_implicit_update ymom_update = domain.quantities['ymomentum'].semi_implicit_update xmom_update = xmom.semi_implicit_update ymom_update = ymom.semi_implicit_update N = len(domain) 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 /= h[k]**(7.0/3) #Update momentum xmom_update[k] += S*uh[k] ymom_update[k] += S*vh[k] def manning_friction_implicit_c(domain): """Wrapper for c version """ #print 'Implicit friction' from shallow_water_ext import manning_friction as manning_friction_c manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update) def manning_friction_explicit(domain): """Apply (Manning) friction to water momentum Wrapper for c version """ # print 'Explicit friction' xmom = domain.quantities['xmomentum'] eta = domain.quantities['friction'].centroid_values xmom_update = xmom.semi_implicit_update ymom_update = ymom.semi_implicit_update xmom_update = xmom.explicit_update ymom_update = ymom.explicit_update N = len(domain) g = domain.g from shallow_water_ext import manning_friction manning_friction(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update) def manning_friction_explicit_c(domain): """Wrapper for c version """ #print 'Explicit friction' xmom = domain.quantities['xmomentum'] ymom = domain.quantities['ymomentum'] w = domain.quantities['stage'].centroid_values z = domain.quantities['elevation'].centroid_values uh = xmom.centroid_values vh = ymom.centroid_values eta = domain.quantities['friction'].centroid_values xmom_update = xmom.explicit_update ymom_update = ymom.explicit_update N = len(domain) eps = domain.minimum_allowed_height g = domain.g from shallow_water_ext import manning_friction manning_friction(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update) from shallow_water_ext import manning_friction as manning_friction_c manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update) # FIXME (Ole): This was implemented for use with one of the analytical solutions (Sampson?) # Is it still needed (30 Oct 2007)? def linear_friction(domain): """Apply linear friction to water momentum #--------------------------------- # Experimental auxiliary functions #--------------------------------- def check_forcefield(f): """Check that f is either except Exception, e: msg = 'Function %s could not be executed:\n%s' %(f, e) #FIXME: Reconsider this semantics # FIXME: Reconsider this semantics raise msg raise msg #Is this really what we want? # Is this really what we want? msg = 'Return vector from function %s ' %f msg += 'must have same lenght as input vectors' phi = args[1] elif len(args) == 1: #Assume vector function returning (s, phi)(t,x,y) # 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] else: #Assume info is in 2 keyword arguments # Assume info is in 2 keyword arguments if len(kwargs) == 2: s_vec = self.speed(t, xc[:,0], xc[:,1]) else: #Assume s is a scalar # Assume s is a scalar try: phi_vec = self.phi(t, xc[:,0], xc[:,1]) else: #Assume phi is a scalar # Assume phi is a scalar try: phi = phi_vec[k] #Convert to radians # Convert to radians phi = phi*pi/180 #Compute velocity vector (u, v) # Compute velocity vector (u, v) u = s*cos(phi) v = s*sin(phi) #Compute wind stress # Compute wind stress S = const * sqrt(u**2 + v**2) xmom_update[k] += S*u #------------------ from anuga.utilities import compile if compile.can_use_C_extension('shallow_water_ext.c'): # Replace python version with c implementations # Underlying C implementations can be accessed from shallow_water_ext import rotate, assign_windfield_values compute_fluxes = compute_fluxes_c extrapolate_second_order_sw=extrapolate_second_order_sw_c gravity = gravity_c manning_friction = manning_friction_implicit_c h_limiter = h_limiter_c balance_deep_and_shallow = balance_deep_and_shallow_c protect_against_infinitesimal_and_negative_heights =\ protect_against_infinitesimal_and_negative_heights_c #distribute_to_vertices_and_edges =\ #              distribute_to_vertices_and_edges_c #(like MH's) else: msg = 'C implementations could not be accessed by %s.\n ' %__file__ msg += 'Make sure compile_all.py has been run as described in ' msg += 'the ANUGA installation guide.' raise Exception, msg
• ## anuga_core/source/anuga/shallow_water/shallow_water_ext.c

 r4733 // Innermost flux function (using stage w=z+h) int flux_function_central(double *q_left, double *q_right, double z_left, double z_right, double n1, double n2, double epsilon, double H0, double g, double *edgeflux, double *max_speed) { int _flux_function_central(double *q_left, double *q_right, double z_left, double z_right, double n1, double n2, double epsilon, double H0, double g, double *edgeflux, double *max_speed) { /*Compute fluxes between volumes for the shallow water wave equation /////////////////////////////////////////////////////////////////// // Gateways to Python PyObject *flux_function_central(PyObject *self, PyObject *args) { // // Gateway to innermost flux function. // This is only used by the unit tests as the c implementation is // normally called by compute_fluxes in this module. PyArrayObject *normal, *ql, *qr,  *edgeflux; double g, epsilon, max_speed, H0, zl, zr; if (!PyArg_ParseTuple(args, "OOOddOddd", &normal, &ql, &qr, &zl, &zr, &edgeflux, &epsilon, &g, &H0)) { PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: flux_function_central could not parse input arguments"); return NULL; } _flux_function_central((double*) ql -> data, (double*) qr -> data, zl, zr, ((double*) normal -> data)[0], ((double*) normal -> data)[1], epsilon, H0, g, (double*) edgeflux -> data, &max_speed); return Py_BuildValue("d", max_speed); } PyObject *gravity(PyObject *self, PyObject *args) { } //Input checks (convert sequences into numeric arrays) // Input checks (convert sequences into numeric arrays) q = (PyArrayObject *) PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0); } //Allocate space for return vector r (don't DECREF) // Allocate space for return vector r (don't DECREF) dimensions[0] = 3; r = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); //Copy // Copy for (i=0; i<3; i++) { ((double *) (r -> data))[i] = ((double *) (q -> data))[i]; } //Get normal and direction // Get normal and direction n1 = ((double *) normal -> data)[0]; n2 = ((double *) normal -> data)[1]; if (direction == -1) n2 = -n2; //Rotate // Rotate _rotate((double *) r -> data, n1, n2); //Release numeric arrays // Release numeric arrays Py_DECREF(q); Py_DECREF(normal); //return result using PyArray to avoid memory leak // Return result using PyArray to avoid memory leak return PyArray_Return(r); } // Edge flux computation (triangle k, edge i) flux_function_central(ql, qr, zl, zr, normals[ki2], normals[ki2+1], epsilon, H0, g, edgeflux, &max_speed); _flux_function_central(ql, qr, zl, zr, normals[ki2], normals[ki2+1], epsilon, H0, g, edgeflux, &max_speed); *xmom_explicit_update, *ymom_explicit_update, *already_computed_flux;//tracks whether the flux across an edge has already been computed //Local variables *already_computed_flux; // Tracks whether the flux across an edge has already been computed // Local variables double timestep, max_speed, epsilon, g, H0; double normal[2], ql[3], qr[3], zl, zr; PyObject *balance_deep_and_shallow(PyObject *self, PyObject *args) { // 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. // //    balance_deep_and_shallow(beta_h, wc, zc, wv, zv, int k, i, n, N, k3; int dimensions[2]; double beta_h; //Safety factor (see config.py) double beta_h; // Safety factor (see config.py) double *hmin, *hmax, hn; n = ((long*) neighbours -> data)[k3+i]; //Initialise hvbar with values from hv // Initialise hvbar with values from hv ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i]; if (n >= 0) { hn = ((double*) hc -> data)[n]; //Neighbour's centroid value hn = ((double*) hc -> data)[n]; // Neighbour's centroid value hmin[k] = min(hmin[k], hn); _limit(N, beta_h, (double*) hc -> data, (double*) hvbar -> data, hmin, hmax); // // //Py_DECREF(domain); //FIXME: NEcessary? // // //Py_DECREF(domain); //FIXME: Necessary? free(hmin); free(hmax); //return result using PyArray to avoid memory leak // Return result using PyArray to avoid memory leak return PyArray_Return(hvbar); //return Py_BuildValue(""); } ////////////////////////////////////////// //------------------------------- // Method table for python module //------------------------------- static struct PyMethodDef MethodTable[] = { /* The cast of the function is necessary since PyCFunction values {"gravity", gravity, METH_VARARGS, "Print out"}, {"manning_friction", manning_friction, METH_VARARGS, "Print out"}, {"flux_function_central", flux_function_central, METH_VARARGS, "Print out"}, {"balance_deep_and_shallow", balance_deep_and_shallow, METH_VARARGS, "Print out"}, {"assign_windfield_values", assign_windfield_values, METH_VARARGS | METH_KEYWORDS, "Print out"}, //{"distribute_to_vertices_and_edges", // distribute_to_vertices_and_edges, METH_VARARGS}, //{"update_conserved_quantities", // update_conserved_quantities, METH_VARARGS}, //{"set_initialcondition", // set_initialcondition, METH_VARARGS}, {NULL, NULL} };
• ## anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

 r4733 from shallow_water_domain import * from shallow_water_domain import flux_function_central as flux_function # Get gateway to C implementation of flux function for direct testing from shallow_water_ext import flux_function_central as flux_function class Weir: #FIXME (Ole): Individual flux tests do NOT test C implementation directly. # Individual flux tests def test_flux_zero_case(self): ql = zeros( 3, Float ) qr = zeros( 3, Float ) normal = zeros( 2, Float ) edgeflux = zeros( 3, Float ) zl = zr = 0. flux, max_speed = flux_function(normal, ql, qr, zl, zr) assert allclose(flux, [0,0,0]) H0 = 0.0 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0) assert allclose(edgeflux, [0,0,0]) assert max_speed == 0. ql = array([w, 0, 0]) qr = array([w, 0, 0]) edgeflux = zeros(3, Float) zl = zr = 0. h = w - (zl+zr)/2 flux, max_speed = flux_function(normal, ql, qr, zl, zr) assert allclose(flux, [0., 0.5*g*h**2, 0.]) H0 = 0.0 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0) assert allclose(edgeflux, [0., 0.5*g*h**2, 0.]) assert max_speed == sqrt(g*h) qr = array([-0.2, 2, 3]) zl = zr = -0.5 flux, max_speed = flux_function(normal, ql, qr, zl, zr) assert allclose(flux, [2.,13.77433333, 20.]) edgeflux = zeros(3, Float) H0 = 0.0 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0) assert allclose(edgeflux, [2.,13.77433333, 20.]) assert allclose(max_speed, 8.38130948661) qr = array([-0.075, 2, 3]) zl = zr = -0.375 flux, max_speed = flux_function(normal, ql, qr, zl, zr) assert allclose(flux, [-3.,-20.0, -30.441]) edgeflux = zeros(3, Float) H0 = 0.0 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0) assert allclose(edgeflux, [-3.,-20.0, -30.441]) assert allclose(max_speed, 11.7146428199) qr = array([-0.075, 2, 3]) zl = zr = -0.375 flux, max_speed = flux_function(normal, ql, qr, zl, zr) assert allclose(flux, [sqrt(2)/2, 4.40221112, 7.3829019]) edgeflux = zeros(3, Float) H0 = 0.0 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0) assert allclose(edgeflux, [sqrt(2)/2, 4.40221112, 7.3829019]) assert allclose(max_speed, 4.0716654239) qr = array([-0.30683287, 0.1071986, 0.05930515]) zl = zr = -0.375 flux, max_speed = flux_function(normal, ql, qr, zl, zr) assert allclose(flux, [-0.04072676, -0.07096636, -0.01604364]) edgeflux = zeros(3, Float) H0 = 0.0 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0) assert allclose(edgeflux, [-0.04072676, -0.07096636, -0.01604364]) assert allclose(max_speed, 1.31414103233) def test_compute_fluxes0(self): #Do a full triangle and check that fluxes cancel out for #the constant stage case # Do a full triangle and check that fluxes cancel out for # the constant stage case a = [0.0, 0.0] assert allclose(domain.neighbour_edges, [[-1,2,-1], [2,0,1], [-1,-1,0],[1,-1,-1]]) zl=zr=0. #Assume flat bed #Flux across right edge of volume 1 zl=zr=0. # Assume flat bed edgeflux = zeros(3, Float) edgeflux0 = zeros(3, Float) edgeflux1 = zeros(3, Float) edgeflux2 = zeros(3, Float) H0 = 0.0 # Flux across right edge of volume 1 normal = domain.get_normal(1,0) ql = domain.get_conserved_quantities(vol_id=1, edge=0) qr = domain.get_conserved_quantities(vol_id=2, edge=2) flux0, max_speed = flux_function(normal, ql, qr, zl, zr) #Check that flux seen from other triangles is inverse max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0, epsilon, g, H0) # Check that flux seen from other triangles is inverse tmp = qr; qr=ql; ql=tmp normal = domain.get_normal(2,2) flux, max_speed = flux_function(normal, ql, qr, zl, zr) assert allclose(flux + flux0, 0.) #Flux across upper edge of volume 1 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0) assert allclose(edgeflux0 + edgeflux, 0.) # Flux across upper edge of volume 1 normal = domain.get_normal(1,1) ql = domain.get_conserved_quantities(vol_id=1, edge=1) qr = domain.get_conserved_quantities(vol_id=3, edge=0) flux1, max_speed = flux_function(normal, ql, qr, zl, zr) #Check that flux seen from other triangles is inverse max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1, epsilon, g, H0) # Check that flux seen from other triangles is inverse tmp = qr; qr=ql; ql=tmp normal = domain.get_normal(3,0) flux, max_speed = flux_function(normal, ql, qr, zl, zr) assert allclose(flux + flux1, 0.) #Flux across lower left hypotenuse of volume 1 max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0) assert allclose(edgeflux1 + edgeflux, 0.) # Flux across lower left hypotenuse of volume 1 normal = domain.get_normal(1,2) ql = domain.get_conserved_quantities(vol_id=1, edge=2) qr = domain.get_conserved_quantities(vol_id=0, edge=1) flux2, max_speed = flux_function(normal, ql, qr, zl, zr) #Check that flux seen from other triangles is inverse max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2, epsilon, g, H0) # Check that flux seen from other triangles is inverse tmp = qr; qr=ql; ql=tmp normal = domain.get_normal(0,1) flux, max_speed = flux_function(normal, ql, qr, zl, zr) assert allclose(flux + flux2, 0.) #Scale by edgelengths, add up anc check that total flux is zero max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0) assert allclose(edgeflux2 + edgeflux, 0.) # Scale by edgelengths, add up anc check that total flux is zero e0 = domain.edgelengths[1, 0] e1 = domain.edgelengths[1, 1] e2 = domain.edgelengths[1, 2] assert allclose(e0*flux0+e1*flux1+e2*flux2, 0.) #Now check that compute_flux yields zeros as well assert allclose(e0*edgeflux0+e1*edgeflux1+e2*edgeflux2, 0.) # Now check that compute_flux yields zeros as well domain.compute_fluxes() zl=zr=0. #Assume flat bed #Flux across right edge of volume 1 edgeflux = zeros(3, Float) edgeflux0 = zeros(3, Float) edgeflux1 = zeros(3, Float) edgeflux2 = zeros(3, Float) H0 = 0.0 # Flux across right edge of volume 1 normal = domain.get_normal(1,0) #Get normal 0 of triangle 1 assert allclose(normal, [1, 0]) qr = domain.get_conserved_quantities(vol_id=2, edge=2) assert allclose(qr, [val2, 0, 0]) flux0, max_speed = flux_function(normal, ql, qr, zl, zr) #Flux across edge in the east direction (as per normal vector) assert allclose(flux0, [-15.3598804, 253.71111111, 0.]) max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0, epsilon, g, H0) # Flux across edge in the east direction (as per normal vector) assert allclose(edgeflux0, [-15.3598804, 253.71111111, 0.]) assert allclose(max_speed, 9.21592824046) normal_opposite = domain.get_normal(2,2) #Get normal 2 of triangle 2 assert allclose(normal_opposite, [-1, 0]) flux_opposite, max_speed = flux_function([-1, 0], ql, qr, zl, zr) assert allclose(flux_opposite, [-15.3598804, -253.71111111, 0.]) max_speed = flux_function(normal_opposite, ql, qr, zl, zr, edgeflux, epsilon, g, H0) #flux_opposite, max_speed = flux_function([-1, 0], ql, qr, zl, zr) assert allclose(edgeflux, [-15.3598804, -253.71111111, 0.]) ql = domain.get_conserved_quantities(vol_id=1, edge=1) qr = domain.get_conserved_quantities(vol_id=3, edge=0) flux1, max_speed = flux_function(normal, ql, qr, zl, zr) assert allclose(flux1, [2.4098563, 0., 123.04444444]) max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1, epsilon, g, H0) assert allclose(edgeflux1, [2.4098563, 0., 123.04444444]) assert allclose(max_speed, 7.22956891292) ql = domain.get_conserved_quantities(vol_id=1, edge=2) qr = domain.get_conserved_quantities(vol_id=0, edge=1) flux2, max_speed = flux_function(normal, ql, qr, zl, zr) assert allclose(flux2, [9.63942522, -61.59685738, -61.59685738]) max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2, epsilon, g, H0) assert allclose(edgeflux2, [9.63942522, -61.59685738, -61.59685738]) assert allclose(max_speed, 7.22956891292) e2 = domain.edgelengths[1, 2] total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1] total_flux = -(e0*edgeflux0+e1*edgeflux1+e2*edgeflux2)/domain.areas[1] assert allclose(total_flux, [-0.68218178, -166.6, -35.93333333]) zl=zr=0 #Assume flat zero bed edgeflux = zeros(3, Float) edgeflux0 = zeros(3, Float) edgeflux1 = zeros(3, Float) edgeflux2 = zeros(3, Float) H0 = 0.0 domain.set_quantity('elevation', zl*ones( (4,3) )) ql = domain.get_conserved_quantities(vol_id=1, edge=0) qr = domain.get_conserved_quantities(vol_id=2, edge=2) flux0, max_speed = flux_function(normal, ql, qr, zl, zr) max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0, epsilon, g, H0) #Flux across upper edge of volume 1 ql = domain.get_conserved_quantities(vol_id=1, edge=1) qr = domain.get_conserved_quantities(vol_id=3, edge=0) flux1, max_speed = flux_function(normal, ql, qr, zl, zr) max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1, epsilon, g, H0) #Flux across lower left hypotenuse of volume 1 ql = domain.get_conserved_quantities(vol_id=1, edge=2) qr = domain.get_conserved_quantities(vol_id=0, edge=1) flux2, max_speed = flux_function(normal, ql, qr, zl, zr) max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2, epsilon, g, H0) #Scale, add up and check that compute_fluxes is correct for vol 1 e2 = domain.edgelengths[1, 2] total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1] total_flux = -(e0*edgeflux0+e1*edgeflux1+e2*edgeflux2)/domain.areas[1] edgeflux = zeros(3, Float) edgeflux0 = zeros(3, Float) edgeflux1 = zeros(3, Float) edgeflux2 = zeros(3, Float) H0 = 0.0 domain.set_quantity('stage', [[val0, val0-1, val0-2], [val1, val1+1, val1], ql = domain.get_conserved_quantities(vol_id=1, edge=0) qr = domain.get_conserved_quantities(vol_id=2, edge=2) flux0, max_speed = flux_function(normal, ql, qr, zl, zr) max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0, epsilon, g, H0) #Flux across upper edge of volume 1 ql = domain.get_conserved_quantities(vol_id=1, edge=1) qr = domain.get_conserved_quantities(vol_id=3, edge=0) flux1, max_speed = flux_function(normal, ql, qr, zl, zr) max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1, epsilon, g, H0) #Flux across lower left hypotenuse of volume 1 ql = domain.get_conserved_quantities(vol_id=1, edge=2) qr = domain.get_conserved_quantities(vol_id=0, edge=1) flux2, max_speed = flux_function(normal, ql, qr, zl, zr) max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2, epsilon, g, H0) #Scale, add up and check that compute_fluxes is correct for vol 1 e2 = domain.edgelengths[1, 2] total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1] total_flux = -(e0*edgeflux0+e1*edgeflux1+e2*edgeflux2)/domain.areas[1] domain.compute_fluxes() if __name__ == "__main__": suite = unittest.makeSuite(Test_Shallow_Water,'test') suite = unittest.makeSuite(Test_Shallow_Water,'test') #suite = unittest.makeSuite(Test_Shallow_Water,'test_extrema') #suite = unittest.makeSuite(Test_Shallow_Water,'test_tight_slope_limiters')
• ## anuga_core/source/anuga/utilities/polygon.py

 r4768 http://www.alienryderflex.com/polygon/ """ Uses underlying C-implementation in polygon_ext.c """ if verbose: print 'Checking input to separate_points_by_polygon' #Input checks assert isinstance(closed, bool), 'Keyword argument "closed" must be boolean' assert isinstance(verbose, bool), 'Keyword argument "verbose" must be boolean' raise msg #if verbose: print 'Checking input to separate_points_by_polygon 2' try: polygon = ensure_numeric(polygon, Float) raise msg #if verbose: print 'check' assert len(polygon.shape) == 2,\ 'Polygon array must be a 2d array of vertices' M = points.shape[0]  #Number of points px = polygon[:,0] py = polygon[:,1] #Used for an optimisation when points are far away from polygon minpx = min(px); maxpx = max(px) minpy = min(py); maxpy = max(py) #Begin main loop indices = zeros(M, Int) inside_index = 0    #Keep track of points inside outside_index = M-1 #Keep track of points outside (starting from end) if verbose: print 'Separating %d points' %M for k in range(M): if verbose: if k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M) #for each point x = points[k, 0] y = points[k, 1] inside = False if not x > maxpx or x < minpx or y > maxpy or y < minpy: #Check polygon for i in range(N): j = (i+1)%N #Check for case where point is contained in line segment if point_on_line(x, y, px[i], py[i], px[j], py[j]): if closed: inside = True else: inside = False break else: #Check if truly inside polygon if py[i] < y and py[j] >= y or\ py[j] < y and py[i] >= y: if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x: inside = not inside if inside: indices[inside_index] = k inside_index += 1 else: indices[outside_index] = k outside_index -= 1 return indices, inside_index def separate_points_by_polygon_c(points, polygon, closed = True, verbose = False): """Determine whether points are inside or outside a polygon C-wrapper """ if verbose: print 'Checking input to separate_points_by_polygon' #Input checks assert isinstance(closed, bool), 'Keyword argument "closed" must be boolean' assert isinstance(verbose, bool), 'Keyword argument "verbose" must be boolean' try: points = ensure_numeric(points, Float) except NameError, e: raise NameError, e except: msg = 'Points could not be converted to Numeric array' raise msg #if verbose: print 'Checking input to separate_points_by_polygon 2' try: polygon = ensure_numeric(polygon, Float) except NameError, e: raise NameError, e except: msg = 'Polygon could not be converted to Numeric array' raise msg #if verbose: print 'check' assert len(polygon.shape) == 2,\ 'Polygon array must be a 2d array of vertices' assert polygon.shape[1] == 2,\ 'Polygon array must have two columns' assert len(points.shape) == 2,\ 'Points array must be a 2d array' assert points.shape[1] == 2,\ 'Points array must have two columns' N = polygon.shape[0] #Number of vertices in polygon M = points.shape[0]  #Number of points from polygon_ext import separate_points_by_polygon if verbose: print 'Allocating array for indices' #if verbose: print 'Calling C-version of inside poly' count = separate_points_by_polygon(points, polygon, indices, int(closed), int(verbose)) from polygon_ext import separate_points_by_polygon as sep_points count = sep_points(points, polygon, indices, int(closed), int(verbose)) if verbose: print 'Found %d points (out of %d) inside polygon'\ from anuga.utilities.compile import can_use_C_extension if can_use_C_extension('polygon_ext.c'): # Underlying C implementations can be accessed from polygon_ext import point_on_line separate_points_by_polygon = separate_points_by_polygon_c else: msg = 'C implementations could not be accessed by %s.\n ' %__file__ msg += 'Make sure compile_all.py has been run as described in ' msg += 'the ANUGA installation guide.' raise Exception, msg
• ## anuga_core/source/anuga/utilities/sparse.py

 r3832 def csr_mv(self, B): """Python version of sparse (CSR) matrix multiplication """ from Numeric import zeros, Float #Assume numeric types from now on if len(B.shape) == 0: #Scalar - use __rmul__ method R = B*self elif len(B.shape) == 1: #Vector assert B.shape[0] == self.N, 'Mismatching dimensions' R = zeros(self.M, Float) #Result #Multiply nonzero elements for i in range(self.M): for ckey in range(self.row_ptr[i],self.row_ptr[i+1]): j = self.colind[ckey] R[i] += self.data[ckey]*B[j] elif len(B.shape) == 2: R = zeros((self.M, B.shape[1]), Float) #Result matrix #Multiply nonzero elements for col in range(R.shape[1]): #For each column for i in range(self.M): for ckey in range(self.row_ptr[i],self.row_ptr[i+1]): j = self.colind[ckey] R[i, col] += self.data[ckey]*B[j,col] else: raise ValueError, 'Dimension too high: d=%d' %len(B.shape) return R #Setup for C extensions # Setup for C extensions from anuga.utilities import compile if compile.can_use_C_extension('sparse_ext.c'): #Replace python version with c implementation # Access underlying c implementations from sparse_ext import csr_mv if __name__ == '__main__': # A little selftest from Numeric import allclose, array, Float
Note: See TracChangeset for help on using the changeset viewer.