Changeset 7823
 Timestamp:
 Jun 11, 2010, 6:14:08 PM (13 years ago)
 Location:
 anuga_work/development/pipeflow
 Files:

 5 added
 5 edited
 2 copied
Legend:
 Unmodified
 Added
 Removed

anuga_work/development/pipeflow/generic_domain.py
r7788 r7823 12 12 13 13 14 class Generic_ Domain:14 class Generic_domain: 15 15 16 16 def __init__(self, … … 587 587 return self.vertices[elem_id,vertex] 588 588 589 def get_area(self, elem_id ):589 def get_area(self, elem_id=None): 590 590 """Return area of element id 591 591 """ 592 592 593 return self.areas[elem_id] 593 if elem_id is None: 594 return sum(self.areas) 595 else: 596 return self.areas[elem_id] 597 594 598 595 599 def get_quantity(self, name, location='vertices', indices = None): 
anuga_work/development/pipeflow/quantity.py
r7784 r7823 1 """Class Quantity  Implements values at each 1d element 1 """ 2 Class Quantity  Implements values at each 1d element 2 3 3 4 To create: … … 17 18 import numpy 18 19 20 21 19 22 class Quantity: 20 23 … … 23 26 #Initialise Quantity using optional vertex values. 24 27 25 from generic_domain import Generic_ Domain28 from generic_domain import Generic_domain 26 29 27 30 msg = 'First argument in Quantity.__init__ ' 28 31 msg += 'must be of class Domain (or a subclass thereof)' 29 assert isinstance(domain, Generic_ Domain), msg32 assert isinstance(domain, Generic_domain), msg 30 33 31 34 if vertex_values is None: … … 52 55 self.centroid_values = numpy.zeros(N, numpy.float) 53 56 self.centroid_backup_values = numpy.zeros(N, numpy.float) 54 #self.edge_values = zeros((N, 2), numpy.float)57 #self.edge_values = numpy.zeros((N, 2), numpy.float) 55 58 #edge values are values of the ends of each interval 56 59 … … 58 61 self.interpolate() 59 62 63 60 64 61 65 #Allocate space for boundary values … … 82 86 """ 83 87 return self.centroid_values.shape[0] 84 88 89 def __neg__(self): 90 """Negate all values in this quantity giving meaning to the 91 expression Q where Q is an instance of class Quantity 92 """ 93 94 Q = Quantity(self.domain) 95 Q.set_values_from_numeric(self.vertex_values) 96 return Q 97 98 99 100 def __add__(self, other): 101 """Add to self anything that could populate a quantity 102 103 E.g other can be a constant, an array, a function, another quantity 104 (except for a filename or points, attributes (for now)) 105  see set_values_from_numeric for details 106 """ 107 108 Q = Quantity(self.domain) 109 Q.set_values_from_numeric(other) 110 111 result = Quantity(self.domain) 112 result.set_values_from_numeric(self.vertex_values + Q.vertex_values) 113 return result 114 115 def __radd__(self, other): 116 """Handle cases like 7+Q, where Q is an instance of class Quantity 117 """ 118 119 return self + other 120 121 def __sub__(self, other): 122 return self + other # Invoke self.__neg__() 123 124 def __mul__(self, other): 125 """Multiply self with anything that could populate a quantity 126 127 E.g other can be a constant, an array, a function, another quantity 128 (except for a filename or points, attributes (for now)) 129  see set_values_from_numeric for details 130 """ 131 132 if isinstance(other, Quantity): 133 Q = other 134 else: 135 Q = Quantity(self.domain) 136 Q.set_values_from_numeric(other) 137 138 result = Quantity(self.domain) 139 140 # The product of vertex_values, edge_values and centroid_values 141 # are calculated and assigned directly without using 142 # set_values_from_numeric (which calls interpolate). Otherwise 143 # centroid values wouldn't be products from q1 and q2 144 result.vertex_values = self.vertex_values * Q.vertex_values 145 result.centroid_values = self.centroid_values * Q.centroid_values 146 147 return result 148 149 def __rmul__(self, other): 150 """Handle cases like 3*Q, where Q is an instance of class Quantity 151 """ 152 153 return self * other 154 155 def __div__(self, other): 156 """Divide self with anything that could populate a quantity 157 158 E.g other can be a constant, an array, a function, another quantity 159 (except for a filename or points, attributes (for now)) 160  see set_values_from_numeric for details 161 162 Zero division is dealt with by adding an epsilon to the divisore 163 FIXME (Ole): Replace this with native INF once we migrate to NumPy 164 """ 165 166 if isinstance(other, Quantity): 167 Q = other 168 else: 169 Q = Quantity(self.domain) 170 Q.set_values_from_numeric(other) 171 172 result = Quantity(self.domain) 173 174 # The quotient of vertex_values, edge_values and centroid_values 175 # are calculated and assigned directly without using 176 # set_values_from_numeric (which calls interpolate). Otherwise 177 # centroid values wouldn't be quotient of q1 and q2 178 result.vertex_values = self.vertex_values/(Q.vertex_values + epsilon) 179 result.centroid_values = self.centroid_values/(Q.centroid_values + epsilon) 180 181 return result 182 183 def __rdiv__(self, other): 184 """Handle cases like 3/Q, where Q is an instance of class Quantity 185 """ 186 187 return self / other 188 189 def __pow__(self, other): 190 """Raise quantity to (numerical) power 191 192 As with __mul__ vertex values are processed entry by entry 193 while centroid and edge values are reinterpolated. 194 195 Example using __pow__: 196 Q = (Q1**2 + Q2**2)**0.5 197 """ 198 199 if isinstance(other, Quantity): 200 Q = other 201 else: 202 Q = Quantity(self.domain) 203 Q.set_values_from_numeric(other) 204 205 result = Quantity(self.domain) 206 207 # The power of vertex_values, edge_values and centroid_values 208 # are calculated and assigned directly without using 209 # set_values_from_numeric (which calls interpolate). Otherwise 210 # centroid values wouldn't be correct 211 result.vertex_values = self.vertex_values ** other 212 result.centroid_values = self.centroid_values ** other 213 214 return result 215 216 def set_values_from_numeric(self, numeric): 217 218 x = numpy.array([1.0,2.0]) 219 y = [1.0,2.0] 220 221 if type(numeric) == type(y): 222 self.set_values_from_array(numeric) 223 elif type(numeric) == type(x): 224 self.set_values_from_array(numeric) 225 elif callable(numeric): 226 self.set_values_from_function(numeric) 227 elif isinstance(numeric, Quantity): 228 self.set_values_from_quantity(numeric) 229 else: # see if it's coercible to a float (float, int or long, etc) 230 try: 231 numeric = float(numeric) 232 except ValueError: 233 msg = ("Illegal type for variable 'numeric': %s" % type(numeric)) 234 raise Exception(msg) 235 self.set_values_from_constant(numeric) 236 237 def set_values_from_constant(self,numeric): 238 239 self.vertex_values[:,:] = numeric 240 self.centroid_values[:,] = numeric 241 242 243 def set_values_from_array(self,numeric): 244 245 self.vertex_values[:,:] = numeric 246 self.interpolate() 247 248 249 def set_values_from_quantity(self,quantity): 250 251 self.vertex_values[:,:] = quantity.vertex_values 252 self.centroid_values[:,] = quantity.centroid_values 253 254 def set_values_from_function(self,function): 255 256 self.vertex_values[:,:] = map(function, self.domain.vertices) 257 self.centroid_values[:,] = map(function, self.domain.centroids) 258 259 85 260 def interpolate(self): 86 261 """ … … 96 271 self.centroid_values[i] = (v0 + v1)/2.0 97 272 273 274 275 98 276 def set_values(self, X, location='vertices'): 99 277 """Set values for quantity 100 278 101 X: Compatible list, numpyarray (see below), constant or function279 X: Compatible list, Numeric array (see below), constant or function 102 280 location: Where values are to be stored. 103 281 Permissible options are: vertices, centroid … … 105 283 106 284 In case of location == 'centroid' the dimension values must 107 be a list of a numpyal array of length N, N being the number285 be a list of a Numerical array of length N, N being the number 108 286 of elements in the mesh. Otherwise it must be of dimension Nx2 109 287 … … 119 297 """ 120 298 121 if location not in ['vertices', 'centroids' ]:122 msg = 'Invalid location: %s, (possible choices vertices, centroids )' %location123 raise msg299 if location not in ['vertices', 'centroids', 'unique_vertices']: 300 msg = 'Invalid location: %s, (possible choices vertices, centroids, unique_vertices)' %location 301 raise Exception(msg) 124 302 125 303 if X is None: 126 304 msg = 'Given values are None' 127 raise msg305 raise Exception(msg) 128 306 129 307 import types … … 134 312 135 313 elif type(X) in [types.FloatType, types.IntType, types.LongType]: 136 if location == 'centroids': 137 self.centroid_values[:] = X 138 else: 139 self.vertex_values[:] = X 314 315 self.centroid_values[:,] = float(X) 316 self.vertex_values[:,:] = float(X) 317 318 elif isinstance(X, Quantity): 319 self.set_array_values(X.vertex_values, location) 140 320 141 321 else: … … 143 323 self.set_array_values(X, location) 144 324 145 if location == 'vertices' :146 #Intialise centroid and edge values325 if location == 'vertices' or location == 'unique_vertices': 326 #Intialise centroid 147 327 self.interpolate() 328 329 if location == 'centroid': 330 self.extrapolate_first_order() 148 331 149 332 … … 176 359 """Set values for quantity 177 360 178 values: numpyarray361 values: Numeric array 179 362 location: Where values are to be stored. 180 Permissible options are: vertices, centroid, edges363 Permissible options are: vertices, centroid, unique_vertices 181 364 Default is "vertices" 182 365 183 366 In case of location == 'centroid' the dimension values must 184 be a list of a numpyal array of length N, N being the number 185 of elements in the mesh. Otherwise it must be of dimension Nx2 367 be a list of a Numerical array of length N, N being the number 368 of elements in the mesh. If location == 'unique_vertices' the 369 dimension values must be a list or a Numeric array of length N+1. 370 Otherwise it must be of dimension Nx2 186 371 187 372 The values will be stored in elements following their … … 195 380 196 381 197 values = numpy.array(values ,numpy.float)382 values = numpy.array(values).astype(numpy.float) 198 383 199 384 N = self.centroid_values.shape[0] 200 385 201 msg = 'Number of values must match number of elements'202 assert values.shape[0] == N, msg203 386 204 387 if location == 'centroids': 388 msg = 'Number of values must match number of elements' 389 assert values.shape[0] == N, msg 205 390 assert len(values.shape) == 1, 'Values array must be 1d' 206 self.centroid_values = values 207 #elif location == 'edges': 208 # assert len(values.shape) == 2, 'Values array must be 2d' 209 # msg = 'Array must be N x 2' 210 # self.edge_values = values 211 else: 391 392 for i in range(N): 393 self.centroid_values[i] = values[i] 394 395 self.vertex_values[:,0] = values 396 self.vertex_values[:,1] = values 397 398 if location == 'vertices': 399 msg = 'Number of values must match number of elements' 400 assert values.shape[0] == N, msg 212 401 assert len(values.shape) == 2, 'Values array must be 2d' 213 402 msg = 'Array must be N x 2' 214 403 assert values.shape[1] == 2, msg 215 404 216 self.vertex_values = values 405 self.vertex_values[:,:] = values[:,:] 406 407 if location == 'unique_vertices': 408 msg = 'Number of values must match number of elements +1' 409 assert values.shape[0] == N+1, msg 410 assert len(values.shape) == 1, 'Values array must be 1d' 411 412 self.vertex_values[:,0] = values[0:1] 413 self.vertex_values[:,1] = values[1:N+1] 414 415 416 217 417 218 418 … … 220 420 """get values for quantity 221 421 222 return X, Compatible list, numpyarray (see below)422 return X, Compatible list, Numeric array (see below) 223 423 location: Where values are to be stored. 224 424 Permissible options are: vertices, centroid … … 226 426 227 427 In case of location == 'centroids' the dimension values must 228 be a list of a numpyal array of length N, N being the number428 be a list of a Numerical array of length N, N being the number 229 429 of elements. Otherwise it must be of dimension Nx3 230 430 … … 236 436 237 437 """ 238 from numpy import take 438 239 439 240 440 if location not in ['vertices', 'centroids', 'unique vertices']: 241 441 msg = 'Invalid location: %s' %location 242 raise msg243 244 import types 442 raise Exception(msg) 443 444 import types, numpy 245 445 assert type(indices) in [types.ListType, types.NoneType, 246 446 numpy.ArrayType],\ … … 262 462 if cells is None: 263 463 msg = 'Unique vertex not associated with cells' 264 raise msg464 raise Exception(msg) 265 465 266 466 # Go through all cells, vertex pairs … … 274 474 if (indices == None): 275 475 indices = range(len(self)) 276 return numpy.take(self.vertex_values,indices)476 return take(self.vertex_values,indices) 277 477 278 478 … … 284 484 """Return vertex values like an OBJ format 285 485 286 The vertex values are returned as one sequence in the 1D numpy.float array A.486 The vertex values are returned as one sequence in the 1D float array A. 287 487 If requested the coordinates will be returned in 1D arrays X. 288 488 … … 308 508 309 509 """ 510 511 310 512 311 513 … … 399 601 """ 400 602 603 604 401 605 N = self.centroid_values.shape[0] 402 606 … … 407 611 denominator = numpy.ones(N, numpy.float)timestep*self.semi_implicit_update 408 612 409 if numpy.sum(numpy.equal(denominator, 0.0)) > 0.0:613 if sum(numpy.equal(denominator, 0.0)) > 0.0: 410 614 msg = 'Zero division in semi implicit update. Call Stephen :)' 411 raise msg615 raise Exception(msg) 412 616 else: 413 617 #Update conserved_quantities from semi implicit updates … … 420 624 """ 421 625 422 #print 'compute_gradient' 626 627 423 628 424 629 N = self.centroid_values.shape[0] … … 510 715 511 716 #print 'compute_minmod_gradients' 717 from numpy import sign 718 512 719 513 720 def xmin(a,b): … … 516 723 def xmic(t,a,b): 517 724 return xmin(t*xmin(a,b), 0.50*(a+b) ) 725 726 518 727 519 728 N = self.centroid_values.shape[0] … … 702 911 """ 703 912 913 704 914 N = self.domain.number_of_elements 705 915 beta = self.domain.beta … … 729 939 730 940 #Deltas between vertex and centroid values 731 dq = zeros(qv.shape, numpy.float)941 dq = numpy.zeros(qv.shape, numpy.float) 732 942 for i in range(2): 733 943 dq[:,i] = qv[:,i]  qc … … 753 963 754 964 from util import minmod, minmod_kurganov, maxmod, vanleer, vanalbada 965 755 966 limiter = self.domain.limiter 756 967 #print limiter 757 968 758 969 #print 'limit_range' 759 N 970 N = self.domain.number_of_elements 760 971 qc = self.centroid_values 761 972 qv = self.vertex_values 762 C 763 X 973 C = self.domain.centroids 974 X = self.domain.vertices 764 975 beta_p = numpy.zeros(N,numpy.float) 765 976 beta_m = numpy.zeros(N,numpy.float) … … 777 988 beta_x[k] = (qc[k+1]qc[k1])/(C[k+1]C[k1]) 778 989 779 dq = zeros(qv.shape, numpy.float)990 dq = numpy.zeros(qv.shape, numpy.float) 780 991 for i in range(2): 781 992 dq[:,i] =self.domain.vertices[:,i]self.domain.centroids … … 817 1028 818 1029 import sys 1030 819 1031 from util import minmod, minmod_kurganov, maxmod, vanleer 820 1032 … … 841 1053 # Check denominator not zero 842 1054 if (qc[k+1]qc[k]) == 0.0: 843 beta_p[k] = numpy.float(sys.maxint)844 beta_m[k] = numpy.float(sys.maxint)1055 beta_p[k] = float(sys.maxint) 1056 beta_m[k] = float(sys.maxint) 845 1057 else: 846 1058 #STEVE LIMIT … … 849 1061 850 1062 #Deltas between vertex and centroid values 851 dq = zeros(qv.shape, numpy.float)1063 dq = numpy.zeros(qv.shape, numpy.float) 852 1064 for i in range(2): 853 1065 dq[:,i] =self.domain.vertices[:,i]self.domain.centroids … … 885 1097 #backup_centroid_values(self) 886 1098 887 self.centroid_backup_values[: ] = (self.centroid_values).astype('f')1099 self.centroid_backup_values[:,] = (self.centroid_values).astype('f') 888 1100 889 1101 def saxpy_centroid_values(self,a,b): 890 1102 # Call correct module function 891 1103 # (either from this module or Cextension) 892 self.centroid_values[: ] = (a*self.centroid_values + b*self.centroid_backup_values).astype('f')1104 self.centroid_values[:,] = (a*self.centroid_values + b*self.centroid_backup_values).astype('f') 893 1105 894 1106 … … 914 1126 if __name__ == "__main__": 915 1127 #from domain import Domain 916 from shallow_water_domain import Domain 917 from numpy import arange 1128 from sww_domain import Domain 918 1129 919 1130 points1 = [0.0, 1.0, 2.0, 3.0] … … 968 1179 raw_input('press return') 969 1180 970 points2 = arange(10)1181 points2 = numpy.arange(10) 971 1182 D2 = Domain(points2) 972 1183 
anuga_work/development/pipeflow/sww_comp_flux_ext.c
r7783 r7823 10 10 11 11 12 /* double max(double a, double b) { */ 13 /* double z; */ 14 /* z=(a>b)?a:b; */ 15 /* return z;} */ 16 17 /* double min(double a, double b) { */ 18 /* double z; */ 19 /* z=(a<b)?a:b; */ 20 /* return z;} */ 21 22 23 // Function to obtain speed from momentum and depth. 24 // This is used by flux functions 25 // Input parameters uh and h may be modified by this function. 26 double _compute_speed(double *uh, 27 double *h, 28 double epsilon, 29 double h0) { 30 31 double u; 32 33 if (*h < epsilon) { 34 *h = 0.0; //Could have been negative 35 u = 0.0; 36 } else { 37 u = *uh/(*h + h0/ *h); 38 } 39 40 41 // Adjust momentum to be consistent with speed 42 *uh = u * *h; 43 44 return u; 45 } 46 47 48 49 12 50 //Innermost flux function (using w=z+h) 13 51 int _flux_function(double *q_left, double *q_right, 52 double z_left, double z_right, 14 53 double normals, double g, double epsilon, double h0, 15 54 double *edgeflux, double *max_speed) { 16 55 17 56 int i; 18 double flux_left[2], flux_right[2];19 double w_left, h_left, uh_left, z_left, u_left, soundspeed_left;20 double w_right, h_right, uh_right, z_right, u_right, soundspeed_right;21 double z,s_max, s_min, denom;57 double ql[2], qr[2], flux_left[2], flux_right[2]; 58 double z, w_left, h_left, uh_left, soundspeed_left, u_left; 59 double w_right, h_right, uh_right, soundspeed_right, u_right; 60 double s_max, s_min, denom; 22 61 23 62 //printf("h0 = %f \n",h0); 24 w_left = q_left[0]; 25 uh_left = q_left[1]*normals; 26 h_left = q_left[2]; 27 z_left = q_left[3]; 28 u_left = q_left[4]*normals; 29 30 w_right = q_right[0]; 31 uh_right = q_right[1]*normals; 32 h_right = q_right[2]; 33 z_right = q_right[3]; 34 u_right = q_right[4]*normals; 63 ql[0] = q_left[0]; 64 ql[1] = q_left[1]; 65 ql[1] = ql[1]*normals; 66 67 qr[0] = q_right[0]; 68 qr[1] = q_right[1]; 69 qr[1] = qr[1]*normals; 35 70 36 71 z = (z_left+z_right)/2.0; 37 72 73 //w_left = ql[0]; 74 //h_left = w_leftz; 75 //uh_left = ql[1]; 76 77 78 79 // Compute speeds in xdirection 80 w_left = ql[0]; 81 h_left = w_leftz; 82 uh_left = ql[1]; 83 84 u_left = _compute_speed(&uh_left, &h_left, epsilon, h0); 85 86 w_right = qr[0]; 87 h_right = w_rightz; 88 uh_right = qr[1]; 89 90 u_right = _compute_speed(&uh_right, &h_right, epsilon, h0); 91 38 92 soundspeed_left = sqrt(g*h_left); 39 93 soundspeed_right = sqrt(g*h_right); … … 48 102 // Flux formulas 49 103 flux_left[0] = u_left*h_left; 50 flux_left[1] = u_left*u _left*h_left + 0.5*g*h_left*h_left;104 flux_left[1] = u_left*uh_left + 0.5*g*h_left*h_left; 51 105 52 106 flux_right[0] = u_right*h_right; 53 flux_right[1] = u_right*u _right*h_right + 0.5*g*h_right*h_right;107 flux_right[1] = u_right*uh_right + 0.5*g*h_right*h_right; 54 108 55 109 // Flux computation … … 60 114 } else { 61 115 edgeflux[0] = s_max*flux_left[0]  s_min*flux_right[0]; 62 edgeflux[0] += s_max*s_min*(q _right[0]q_left[0]);116 edgeflux[0] += s_max*s_min*(qr[0]ql[0]); 63 117 edgeflux[0] /= denom; 64 118 edgeflux[1] = s_max*flux_left[1]  s_min*flux_right[1]; 65 edgeflux[1] += s_max*s_min*(q _right[1]q_left[1]);119 edgeflux[1] += s_max*s_min*(qr[1]ql[1]); 66 120 edgeflux[1] /= denom; 67 121 edgeflux[1] *= normals; 68 122 69 70 123 // Maximal wavespeed 124 *max_speed = max(fabs(s_max), fabs(s_min)); 71 125 } 72 73 }126 return 0; 127 } 74 128 75 129 … … 90 144 double* xmom_edge_values, 91 145 double* bed_edge_values, 92 double* height_edge_values,93 double* velocity_edge_values,94 146 double* stage_boundary_values, 95 147 double* xmom_boundary_values, 96 double* bed_boundary_values,97 double* height_boundary_values,98 double* velocity_boundary_values,99 148 double* stage_explicit_update, 100 149 double* xmom_explicit_update, … … 102 151 double* max_speed_array) { 103 152 104 double flux[2], ql[ 5], qr[5], edgeflux[2];105 double max_speed, normal;153 double flux[2], ql[2], qr[2], edgeflux[2]; 154 double zl, zr, max_speed, normal; 106 155 int k, i, ki, n, m, nm=0; 107 156 … … 116 165 ql[0] = stage_edge_values[ki]; 117 166 ql[1] = xmom_edge_values[ki]; 118 ql[2] = bed_edge_values[ki]; 119 ql[3] = height_edge_values[ki]; 120 ql[4] = velocity_edge_values[ki]; 167 zl = bed_edge_values[ki]; 121 168 122 169 n = neighbours[ki]; … … 125 172 qr[0] = stage_boundary_values[m]; 126 173 qr[1] = xmom_boundary_values[m]; 127 qr[2] = bed_edge_values[m]; 128 qr[3] = height_edge_values[m]; 129 qr[4] = velocity_edge_values[m]; 174 zr = zl; 130 175 } else { 131 176 m = neighbour_vertices[ki]; … … 133 178 qr[0] = stage_edge_values[nm]; 134 179 qr[1] = xmom_edge_values[nm]; 135 qr[2] = bed_edge_values[nm]; 136 qr[3] = height_edge_values[nm]; 137 qr[4] = velocity_edge_values[nm]; 180 zr = bed_edge_values[nm]; 138 181 } 139 182 140 183 normal = normals[ki]; 141 _flux_function(ql, qr, normal, g, epsilon, h0, edgeflux, &max_speed);184 _flux_function(ql, qr, zl, zr, normal, g, epsilon, h0, edgeflux, &max_speed); 142 185 flux[0] = edgeflux[0]; 143 186 flux[1] = edgeflux[1]; … … 152 195 } 153 196 } 154 197 } // End edge i (and neighbour n) 155 198 flux[0] /= areas[k]; 156 199 stage_explicit_update[k] = flux[0]; … … 163 206 return timestep; 164 207 } 208 209 210 211 212 213 214 165 215 166 216 //========================================================================= … … 168 218 //========================================================================= 169 219 PyObject *compute_fluxes_ext(PyObject *self, PyObject *args) { 170 171 PyObject 220 221 PyArrayObject 222 *neighbours, 223 *neighbour_vertices, 224 *normals, 225 *areas, 226 *stage_edge_values, 227 *xmom_edge_values, 228 *bed_edge_values, 229 *stage_boundary_values, 230 *xmom_boundary_values, 231 *stage_explicit_update, 232 *xmom_explicit_update, 233 *max_speed_array; 234 235 double timestep, epsilon, g, h0, cfl; 236 int number_of_elements; 237 238 // Convert Python arguments to C 239 if (!PyArg_ParseTuple(args, "dddddOOOOOOOOOOOiO", 240 &cfl, 241 ×tep, 242 &epsilon, 243 &g, 244 &h0, 245 &neighbours, 246 &neighbour_vertices, 247 &normals, 248 &areas, 249 &stage_edge_values, 250 &xmom_edge_values, 251 &bed_edge_values, 252 &stage_boundary_values, 253 &xmom_boundary_values, 254 &stage_explicit_update, 255 &xmom_explicit_update, 256 &number_of_elements, 257 &max_speed_array)) { 258 PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext could not parse input"); 259 return NULL; 260 } 261 262 263 // Call underlying flux computation routine and update 264 // the explicit update arrays 265 timestep = _compute_fluxes_ext(cfl, 266 timestep, 267 epsilon, 268 g, 269 h0, 270 (long*) neighbours > data, 271 (long*) neighbour_vertices > data, 272 (double*) normals > data, 273 (double*) areas > data, 274 (double*) stage_edge_values > data, 275 (double*) xmom_edge_values > data, 276 (double*) bed_edge_values > data, 277 (double*) stage_boundary_values > data, 278 (double*) xmom_boundary_values > data, 279 (double*) stage_explicit_update > data, 280 (double*) xmom_explicit_update > data, 281 number_of_elements, 282 (double*) max_speed_array > data); 283 284 // Return updated flux timestep 285 return Py_BuildValue("d", timestep); 286 } 287 288 289 // 290 PyObject *compute_fluxes_ext_short(PyObject *self, PyObject *args) { 291 292 PyObject 172 293 *domain, 173 *stage, 174 *xmom, 175 *bed, 176 *height, 177 *velocity; 178 179 PyArrayObject 180 *neighbours, 294 *stage, 295 *xmom, 296 *bed; 297 298 PyArrayObject 299 *neighbours, 181 300 *neighbour_vertices, 182 *normals, 301 *normals, 183 302 *areas, 184 303 *stage_vertex_values, 185 304 *xmom_vertex_values, 186 305 *bed_vertex_values, 187 *height_vertex_values,188 *velocity_vertex_values,189 306 *stage_boundary_values, 190 307 *xmom_boundary_values, 191 *bed_boundary_values,192 *height_boundary_values,193 *velocity_boundary_values,194 308 *stage_explicit_update, 195 309 *xmom_explicit_update, 196 310 *max_speed_array; 197 311 198 312 double timestep, epsilon, g, h0, cfl; 199 313 int number_of_elements; 200 314 201 315 202 316 // Convert Python arguments to C 203 if (!PyArg_ParseTuple(args, "dOOOO OO",317 if (!PyArg_ParseTuple(args, "dOOOO", 204 318 ×tep, 205 319 &domain, 206 320 &stage, 207 321 &xmom, 208 &bed, 209 &height, 210 &velocity)) { 322 &bed)) { 211 323 PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext_short could not parse input"); 212 324 return NULL; … … 217 329 g = get_python_double(domain,"g"); 218 330 h0 = get_python_double(domain,"h0"); 219 cfl = get_python_double(domain," cfl");220 221 331 cfl = get_python_double(domain,"CFL"); 332 333 222 334 neighbours = get_consecutive_array(domain, "neighbours"); 223 neighbour_vertices= get_consecutive_array(domain, "neighbour_vertices"); 335 neighbour_vertices= get_consecutive_array(domain, "neighbour_vertices"); 224 336 normals = get_consecutive_array(domain, "normals"); 225 areas = get_consecutive_array(domain, "areas"); 337 areas = get_consecutive_array(domain, "areas"); 226 338 max_speed_array = get_consecutive_array(domain, "max_speed_array"); 227 228 stage_vertex_values = get_consecutive_array(stage, "vertex_values"); 229 xmom_vertex_values = get_consecutive_array(xmom, "vertex_values"); 230 bed_vertex_values = get_consecutive_array(bed, "vertex_values"); 231 height_vertex_values = get_consecutive_array(height, "vertex_values"); 232 velocity_vertex_values = get_consecutive_array(velocity, "vertex_values"); 233 234 stage_boundary_values = get_consecutive_array(stage, "boundary_values"); 235 xmom_boundary_values = get_consecutive_array(xmom, "boundary_values"); 236 bed_boundary_values = get_consecutive_array(bed, "boundary_values"); 237 height_boundary_values = get_consecutive_array(height, "boundary_values"); 238 velocity_boundary_values = get_consecutive_array(velocity, "boundary_values"); 239 240 241 stage_explicit_update = get_consecutive_array(stage, "explicit_update"); 242 xmom_explicit_update = get_consecutive_array(xmom, "explicit_update"); 339 340 stage_vertex_values = get_consecutive_array(stage, "vertex_values"); 341 xmom_vertex_values = get_consecutive_array(xmom, "vertex_values"); 342 bed_vertex_values = get_consecutive_array(bed, "vertex_values"); 343 344 stage_boundary_values = get_consecutive_array(stage, "boundary_values"); 345 xmom_boundary_values = get_consecutive_array(xmom, "boundary_values"); 346 347 348 stage_explicit_update = get_consecutive_array(stage, "explicit_update"); 349 xmom_explicit_update = get_consecutive_array(xmom, "explicit_update"); 350 351 243 352 244 353 number_of_elements = stage_vertex_values > dimensions[0]; 245 354 246 // Call underlying flux computation routine and update 247 // the explicit update arrays 355 356 357 // Call underlying flux computation routine and update 358 // the explicit update arrays 248 359 timestep = _compute_fluxes_ext( 249 360 cfl, … … 259 370 (double*) xmom_vertex_values > data, 260 371 (double*) bed_vertex_values > data, 261 (double*) height_vertex_values > data,262 (double*) velocity_vertex_values > data,263 372 (double*) stage_boundary_values > data, 264 373 (double*) xmom_boundary_values > data, 265 (double*) bed_boundary_values > data,266 (double*) height_boundary_values > data,267 (double*) velocity_boundary_values > data,268 374 (double*) stage_explicit_update > data, 269 375 (double*) xmom_explicit_update > data, … … 279 385 Py_DECREF(xmom_vertex_values); 280 386 Py_DECREF(bed_vertex_values); 281 Py_DECREF(height_vertex_values);282 Py_DECREF(velocity_vertex_values);283 387 Py_DECREF(stage_boundary_values); 284 388 Py_DECREF(xmom_boundary_values); 285 Py_DECREF(bed_boundary_values);286 Py_DECREF(height_boundary_values);287 Py_DECREF(velocity_boundary_values);288 389 Py_DECREF(stage_explicit_update); 289 390 Py_DECREF(xmom_explicit_update); … … 291 392 292 393 394 395 293 396 // Return updated flux timestep 294 397 return Py_BuildValue("d", timestep); 295 398 } 296 399 400 297 401 298 402 … … 303 407 static struct PyMethodDef MethodTable[] = { 304 408 {"compute_fluxes_ext", compute_fluxes_ext, METH_VARARGS, "Print out"}, 409 {"compute_fluxes_ext_short", compute_fluxes_ext_short, METH_VARARGS, "Print out"}, 305 410 {NULL, NULL} 306 411 }; 
anuga_work/development/pipeflow/sww_domain.py
r7781 r7823 50 50 51 51 #Shallow water domain 52 class Domain(Generic_ Domain):52 class Domain(Generic_domain): 53 53 54 54 def __init__(self, coordinates, boundary = None, tagged_elements = None): … … 57 57 evolved_quantities = ['stage', 'xmomentum', 'elevation', 'height', 'velocity'] 58 58 other_quantities = ['friction'] 59 Generic_ Domain.__init__(self,59 Generic_domain.__init__(self, 60 60 coordinates = coordinates, 61 61 boundary = boundary, … … 143 143 144 144 145 timestep = numpy.float(sys.maxint)145 timestep = float(sys.maxint) 146 146 147 147 stage = self.quantities['stage'] … … 152 152 153 153 154 from sww_comp_flux_ext import compute_fluxes_ext 155 156 self.flux_timestep = compute_fluxes_ext(timestep,self,stage,xmom,bed,height,velocity) 154 from sww_comp_flux_ext import compute_fluxes_ext_short 155 156 #self.flux_timestep = compute_fluxes_ext(timestep,self,stage,xmom,bed,height,velocity) 157 158 159 self.flux_timestep = compute_fluxes_ext_short(timestep,self,stage,xmom,bed) 160 157 161 158 162 … … 168 172 169 173 #Call basic machinery from parent class 170 for t in Generic_ Domain.evolve(self, yieldstep, finaltime,duration,174 for t in Generic_domain.evolve(self, yieldstep, finaltime,duration, 171 175 skip_initial_step): 172 176 … … 266 270 ## else: 267 271 268 u_C[: ] = uh_C/(h_C + h0/h_C)272 u_C[:,] = uh_C/(h_C + h0/h_C) 269 273 270 274 for name in [ 'velocity', 'stage' ]: … … 286 290 xmom_V = domain.quantities['xmomentum'].vertex_values 287 291 288 h_V[: ] = stage_V  bed_V292 h_V[:,:] = stage_V  bed_V 289 293 for i in range(len(h_C)): 290 294 for j in range(2): … … 297 301 stage_V[i,(j+1)%2] = stage_V[i,(j+1)%2] + dh 298 302 299 xmom_V[: ] = u_V * h_V303 xmom_V[:,:] = u_V * h_V 300 304 301 305 return … … 314 318 315 319 #Shortcuts 316 wc = domain.quantities['stage'].centroid_values317 zc = domain.quantities['elevation'].centroid_values320 wc = domain.quantities['stage'].centroid_values 321 zc = domain.quantities['elevation'].centroid_values 318 322 xmomc = domain.quantities['xmomentum'].centroid_values 319 323 hc = wc  zc #Water depths at centroids 
anuga_work/development/pipeflow/sww_python.py
r7789 r7823 58 58 for k in range(N): 59 59 60 flux[: ] = 0. #Reset work numpy.array60 flux[:,] = 0. #Reset work numpy.array 61 61 #for i in range(3): 62 62 for i in range(2):
Note: See TracChangeset
for help on using the changeset viewer.