Changeset 7842
 Timestamp:
 Jun 15, 2010, 5:34:24 PM (10 years ago)
 Location:
 anuga_work/development/2010projects/anuga_1d
 Files:

 2 added
 5 edited
Legend:
 Unmodified
 Added
 Removed

anuga_work/development/2010projects/anuga_1d/channel/channel_domain.py
r7839 r7842 3 3 the shallow water wave equation. 4 4 5 This module contains a specialisation of class Domain from module domain.py 6 consisting of methods specific to the Shallow Water Wave Equation 5 This module contains a specialisation of class Generic_domain from module 6 generic_domain.py 7 consisting of methods specific to Channel flow using the Shallow Water Wave Equation 7 8 8 9 This particular modification of the Domain class implements the ability to … … 14 15 15 16 where 16 !!!! NOTE THIS NEEDS UPDATING !!!! 17 U = [A, Q] 18 E = [Q, Q^2/A + g h^2/2]17 18 U = [A, Q] = [b*h, u*b*h] 19 E = [Q, Q^2/A + g*b*h^2/2] 19 20 S represents source terms forcing the system 20 (e.g. gravity, friction, wind stress, ...) 21 (e.g. gravity, boundary_stree, friction, wind stress, ...) 22 gravity = g*b*h*z_x 23 boundary_stress = 1/2*g*b_x*h^2 21 24 22 25 and _t, _x, _y denote the derivative with respect to t, x and y respectiely. … … 25 28 26 29 symbol variable name explanation 30 A area Wetted area = b*h 31 Q discharge flux of water = u*b*h 27 32 x x horizontal distance from origin [m] 28 33 z elevation elevation of bed on which flow is modelled [m] … … 31 36 u speed in the x direction [m/s] 32 37 uh xmomentum momentum in the x direction [m^2/s] 33 38 b width width of channel 34 39 eta mannings friction coefficient [to appear] 35 40 nu wind stress coefficient [to appear] 36 41 37 The conserved quantities are w, uh42 The conserved quantities are A, Q 38 43  39 44 For details see e.g. … … 48 53 49 54 50 from domain import * 51 Generic_Domain = Domain #Rename 55 from anuga_1d.generic.generic_domain import * 56 import numpy 57 52 58 53 59 #Shallow water domain 54 class Domain(Generic_ Domain):60 class Domain(Generic_domain): 55 61 56 62 def __init__(self, coordinates, boundary = None, tagged_elements = None): … … 59 65 evolved_quantities = ['area', 'discharge', 'elevation', 'height', 'velocity','width','stage'] 60 66 other_quantities = ['friction'] 61 Generic_ Domain.__init__(self,67 Generic_domain.__init__(self, 62 68 coordinates = coordinates, 63 69 boundary = boundary, … … 67 73 tagged_elements = tagged_elements) 68 74 69 from config import minimum_allowed_height, g, h075 from anuga_1d.config import minimum_allowed_height, g, h0 70 76 self.minimum_allowed_height = minimum_allowed_height 71 77 self.g = g … … 74 80 self.discontinousb = False 75 81 76 82 83 # forcing terms gravity and boundary stress are included in the flux calculation 77 84 #self.forcing_terms.append(gravity) 78 85 #self.forcing_terms.append(boundary_stress) … … 88 95 89 96 #Reduction operation for get_vertex_values 90 from util import mean97 from anuga_1d.generic.util import mean 91 98 self.reduction = mean 92 99 #self.reduction = min #Looks better near steep slopes … … 94 101 self.set_quantities_to_be_stored(['area','discharge']) 95 102 96 self.__doc__ = 'channel_domain _Ab'103 self.__doc__ = 'channel_domain' 97 104 98 105 self.check_integrity() … … 118 125 msg = 'Fifth evolved quantity must be "velocity"' 119 126 assert self.evolved_quantities[4] == 'velocity', msg 120 msg = ' Fifth evolved quantity must be "width"'127 msg = 'Sixth evolved quantity must be "width"' 121 128 assert self.evolved_quantities[5] == 'width', msg 122 123 Generic_Domain.check_integrity(self) 129 msg = 'Seventh evolved quantity must be "stage"' 130 assert self.evolved_quantities[6] == 'stage', msg 131 132 Generic_domain.check_integrity(self) 124 133 125 134 def compute_fluxes(self): … … 140 149 # 141 150 def compute_fluxes_channel(domain): 142 from Numeric import zeros, Float143 151 import sys 144 145 146 152 timestep = float(sys.maxint) 147 153 148 area = domain.quantities['area']149 discharge 150 bed = domain.quantities['elevation']151 height = domain.quantities['height']152 velocity = domain.quantities['velocity']153 width = domain.quantities['width']154 155 156 from channel_domain_ext import compute_fluxes_channel_ext154 area = domain.quantities['area'] 155 discharge = domain.quantities['discharge'] 156 bed = domain.quantities['elevation'] 157 height = domain.quantities['height'] 158 velocity = domain.quantities['velocity'] 159 width = domain.quantities['width'] 160 161 162 from anuga_1d.channel.channel_domain_ext import compute_fluxes_channel_ext 157 163 domain.flux_timestep = compute_fluxes_channel_ext(timestep,domain,area,discharge,bed,height,velocity,width) 158 164 … … 168 174 169 175 import sys 170 from Numeric import zeros, Float 171 from config import epsilon, h0 172 ## linearb(domain) 173 174 176 from anuga_1d.config import epsilon, h0 175 177 176 178 … … 178 180 179 181 #Shortcuts 180 Area= domain.quantities['area']181 Discharge = domain.quantities['discharge']182 Bed= domain.quantities['elevation']183 Height= domain.quantities['height']184 Velocity= domain.quantities['velocity']185 Width= domain.quantities['width']186 Stage= domain.quantities['stage']182 area = domain.quantities['area'] 183 discharge = domain.quantities['discharge'] 184 bed = domain.quantities['elevation'] 185 height = domain.quantities['height'] 186 velocity = domain.quantities['velocity'] 187 width = domain.quantities['width'] 188 stage = domain.quantities['stage'] 187 189 188 190 #Arrays 189 a_C = Area.centroid_values190 d_C = Discharge.centroid_values191 z_C = Bed.centroid_values192 h_C = Height.centroid_values193 u_C = Velocity.centroid_values194 b_C = Width.centroid_values195 w_C = Stage.centroid_values191 a_C = area.centroid_values 192 d_C = discharge.centroid_values 193 z_C = bed.centroid_values 194 h_C = height.centroid_values 195 u_C = velocity.centroid_values 196 b_C = width.centroid_values 197 w_C = stage.centroid_values 196 198 197 199 if domain.setstageflag: 198 for i in range(len(a_C)): 199 a_C[i]=(w_C[i]z_C[i])*b_C[i] 200 200 a_C[:,] = (w_Cz_C)*b_C 201 201 domain.setstageflag = False 202 202 … … 205 205 206 206 207 h0 = 1.0e12 208 #print id(h_C) 209 for i in range(N): 210 211 if a_C[i] <= h0: 212 a_C[i] = 0.0 213 h_C[i] = 0.0 214 d_C[i] = 0.0 215 u_C[i] = 0.0 216 w_C[i] = z_C[i] 217 218 219 220 else: 221 222 if b_C[i]<=h0: 223 a_C[i] = 0.0 224 h_C[i] = 0.0 225 d_C[i] = 0.0 226 u_C[i] = 0.0 227 w_C[i] = z_C[i] 228 229 else: 230 h_C[i] = a_C[i]/(b_C[i]+h0/b_C[i]) 231 w_C[i] = h_C[i]+z_C[i] 232 u_C[i] = d_C[i]/(a_C[i]+h0/a_C[i]) 207 # FIXME (SGR): Replace with numpy.where to speed up code 208 h0 = 1.0e12 209 210 h_C[:] = numpy.where( (a_C>h0)  (b_C>h0), a_C/(b_C+h0/b_C), 0.0 ) 211 u_C[:] = numpy.where( (a_C>h0)  (b_C>h0), d_C/(a_C+h0/a_C), 0.0 ) 212 213 a_C[:] = numpy.where( (a_C>h0)  (b_C>h0), a_C, 0.0 ) 214 b_C[:] = numpy.where( (a_C>h0)  (b_C>h0), b_C, 0.0 ) 215 d_C[:] = numpy.where( (a_C>h0)  (b_C>h0), d_C, 0.0 ) 216 217 w_C[:] = h_C + z_C 218 219 # for i in range(N): 220 # 221 # if a_C[i] <= h0: 222 # a_C[i] = 0.0 223 # h_C[i] = 0.0 224 # d_C[i] = 0.0 225 # u_C[i] = 0.0 226 # w_C[i] = z_C[i] 227 # else: 228 # 229 # if b_C[i]<=h0: 230 # a_C[i] = 0.0 231 # h_C[i] = 0.0 232 # d_C[i] = 0.0 233 # u_C[i] = 0.0 234 # w_C[i] = z_C[i] 235 # 236 # else: 237 # h_C[i] = a_C[i]/(b_C[i]+h0/b_C[i]) 238 # w_C[i] = h_C[i]+z_C[i] 239 # u_C[i] = d_C[i]/(a_C[i]+h0/a_C[i]) 233 240 234 235 236 241 242 bed.extrapolate_second_order() 243 237 244 for name in ['velocity','stage']: 238 245 Q = domain.quantities[name] … … 243 250 else: 244 251 raise 'Unknown order' 245 a_V = domain.quantities['area'].vertex_values 246 w_V = domain.quantities['stage'].vertex_values 247 z_V = domain.quantities['elevation'].vertex_values 248 h_V = domain.quantities['height'].vertex_values 249 u_V = domain.quantities['velocity'].vertex_values 250 d_V = domain.quantities['discharge'].vertex_values 251 b_V = domain.quantities['width'].vertex_values 252 253 254 a_V = area.vertex_values 255 w_V = stage.vertex_values 256 z_V = bed.vertex_values 257 h_V = height.vertex_values 258 u_V = velocity.vertex_values 259 d_V = discharge.vertex_values 260 b_V = width.vertex_values 252 261 253 262 254 255 for i in range(len(h_C)): 256 for j in range(2): 257 ## if b_V[i,j] < h0 : 258 ## a_V[i,j]=0 259 ## h_V[i,j]=0 260 ## d_V[i,j]=0 261 ## u_V[i,j]=0 262 ## else: 263 264 h_V[i,j] = w_V[i,j]z_V[i,j] 265 if h_V[i,j]<h0: 266 h_V[i,j]=0 267 w_V[i,j]=z_V[i,j] 268 a_V[i,j] = b_V[i,j]*h_V[i,j] 269 d_V[i,j]=u_V[i,j]*a_V[i,j] 270 271 272 273 274 263 #FIXME (SGR): replace with numpy.where 264 265 266 h_V[:] = w_Vz_V 267 268 w_V[:] = numpy.where(h_V > h0, w_V, z_V) 269 h_V[:] = numpy.where(h_V > h0, h_V, 0.0) 270 a_V[:] = b_V*h_V 271 d_V[:] = u_V*a_V 272 273 # for i in range(len(h_C)): 274 # for j in range(2): 275 # h_V[i,j] = w_V[i,j]z_V[i,j] 276 # if h_V[i,j]<h0: 277 # h_V[i,j]=0 278 # w_V[i,j]=z_V[i,j] 279 # a_V[i,j] = b_V[i,j]*h_V[i,j] 280 # d_V[i,j]=u_V[i,j]*a_V[i,j] 275 281 276 282 return … … 278 284 279 285 # 280 #Boundaries  specific to the shallow_water_vel_domain286 #Boundaries  specific to the channel_domain 281 287 # 282 288 class Reflective_boundary(Boundary): … … 298 304 #Handy shorthands 299 305 self.normals = domain.normals 300 self.area = domain.quantities['area'].vertex_values306 self.area = domain.quantities['area'].vertex_values 301 307 self.discharge = domain.quantities['discharge'].vertex_values 302 308 self.bed = domain.quantities['elevation'].vertex_values … … 306 312 self.stage = domain.quantities['stage'].vertex_values 307 313 308 from Numeric import zeros, Float 309 #self.conserved_quantities = zeros(3, Float) 310 self.evolved_quantities = zeros(7, Float) 314 self.evolved_quantities = numpy.zeros(7, numpy.float) 311 315 312 316 def __repr__(self): … … 320 324 321 325 q = self.evolved_quantities 322 q[0] = self.area[vol_id, edge_id]326 q[0] = self.area[vol_id, edge_id] 323 327 q[1] = self.discharge[vol_id, edge_id] 324 q[2] = self.bed[vol_id, edge_id]325 q[3] = self.height[vol_id, edge_id]328 q[2] = self.bed[vol_id, edge_id] 329 q[3] = self.height[vol_id, edge_id] 326 330 q[4] = self.velocity[vol_id, edge_id] 327 q[5] = self.width[vol_id,edge_id] 328 q[6] = self.stage[vol_id,edge_id] 329 330 #print "In Reflective q ",q 331 331 q[5] = self.width[vol_id,edge_id] 332 q[6] = self.stage[vol_id,edge_id] 332 333 333 334 return q … … 348 349 raise msg 349 350 350 from Numeric import array, Float 351 self.evolved_quantities=array(evolved_quantities).astype(Float) 351 assert len(evolved_quantities) == 7 352 353 self.evolved_quantities=numpy.array(evolved_quantities,numpy.float) 352 354 353 355 def __repr__(self): … … 370 372 371 373 372 Area = domain.quantities['area']374 Area = domain.quantities['area'] 373 375 Discharge = domain.quantities['discharge'] 374 376 Elevation = domain.quantities['elevation'] … … 402 404 def boundary_stress(domain): 403 405 404 405 406 from util import gradient 406 407 from Numeric import zeros, Float, array, sum 407 408 409 408 410 409 Area = domain.quantities['area'] … … 415 414 416 415 discharge_ud = Discharge.explicit_update 417 418 419 416 420 417 h = Height.vertex_values 421 418 b = Width.vertex_values … … 451 448 452 449 uh = domain.quantities['xmomentum'].centroid_values 453 #vh = domain.quantities['ymomentum'].centroid_values454 450 eta = domain.quantities['friction'].centroid_values 455 451 456 452 xmom_update = domain.quantities['xmomentum'].semi_implicit_update 457 #ymom_update = domain.quantities['ymomentum'].semi_implicit_update458 453 459 454 N = domain.number_of_elements … … 502 497 503 498 504 def check_forcefield(f):505 """Check that f is either506 1: a callable object f(t,x,y), where x and y are vectors507 and that it returns an array or a list of same length508 as x and y509 2: a scalar510 """511 512 from Numeric import ones, Float, array513 514 515 if callable(f):516 #N = 3517 N = 2518 #x = ones(3, Float)519 #y = ones(3, Float)520 x = ones(2, Float)521 #y = ones(2, Float)522 523 try:524 #q = f(1.0, x=x, y=y)525 q = f(1.0, x=x)526 except Exception, e:527 msg = 'Function %s could not be executed:\n%s' %(f, e)528 #FIXME: Reconsider this semantics529 raise msg530 531 try:532 q = array(q).astype(Float)533 except:534 msg = 'Return value from vector function %s could ' %f535 msg += 'not be converted into a Numeric array of floats.\n'536 msg += 'Specified function should return either list or array.'537 raise msg538 539 #Is this really what we want?540 msg = 'Return vector from function %s ' %f541 msg += 'must have same lenght as input vectors'542 assert len(q) == N, msg543 544 else:545 try:546 f = float(f)547 except:548 msg = 'Force field %s must be either a scalar' %f549 msg += ' or a vector function'550 raise msg551 return f552 499 553 500 def linearb(domain): 
anuga_work/development/2010projects/anuga_1d/channel/channel_domain_ext.c
r7839 r7842 1 1 #include "Python.h" 2 #include " Numeric/arrayobject.h"2 #include "numpy/arrayobject.h" 3 3 #include "math.h" 4 4 #include <stdio.h> 
anuga_work/development/2010projects/anuga_1d/channel/channel_flow_1_padarn.py
r7839 r7842 2 2 from math import sqrt, pow, pi 3 3 4 from Numeric import allclose, array, zeros, ones, Float, take, sqrt 4 import numpy 5 5 6 from channel_domain import * 7 from config import g, epsilon 6 from anuga_1d.channel.channel_domain import * 7 from anuga_1d.config import g, epsilon 8 from anuga_1d.generic.generic_mesh import interval_mesh 8 9 9 10 … … 21 22 22 23 def bed(x): 23 y = zeros(len(x),Float) 24 for i in range(len(x)): 25 if x[i]<525 and x[i]>475: 26 y[i]=1 27 else: 28 y[i]=0 29 return y 24 y = numpy.ones(len(x),numpy.float) 25 26 return numpy.where( (x<525) & (x>475),y,0.0) 27 30 28 31 29 def initial_discharge(x): … … 35 33 36 34 # Set final time and yield time for simulation 37 finaltime = 10.038 yieldstep = finaltime35 finaltime = 50.0 36 yieldstep = 10.0 39 37 40 38 # Length of channel (m) 41 39 L = 1000.0 42 40 # Define the number of cells 43 number_of_cells = [200] 44 45 # Define cells for finite volume and their size 46 N = int(number_of_cells[0]) 47 print "Evaluating domain with %d cells" %N 48 cell_len = L/N # Origin = 0.0 49 points = zeros(N+1,Float) 50 51 # Define the centroid points 52 for j in range(N+1): 53 points[j] = j*cell_len 41 N = 200 54 42 55 43 # Create domain with centroid points as defined above 56 domain = Domain(points) 44 domain = Domain(*interval_mesh(N)) 45 57 46 58 47 # Set initial values of quantities  default to zero … … 63 52 64 53 # Set boundry type, order, timestepping method and limiter 65 domain.set_boundary({'exterior': Dirichlet_boundary([14,20,0,1.4,20/14,9])}) 54 Bd = Dirichlet_boundary([14,20,0,1.4,20/14,9,1.4]) 55 domain.set_boundary({'left': Bd , 'right' : Bd }) 56 57 66 58 domain.order = 2 67 59 domain.set_timestepping_method('rk2') … … 76 68 domain.write_time() 77 69 78 N = float(N) 79 HeightC = domain.quantities['height'].centroid_values 70 71 exit() 72 73 HeightC = domain.quantities['height'].centroid_values 80 74 DischargeC = domain.quantities['discharge'].centroid_values 81 C = domain.centroids75 C = domain.centroids 82 76 print 'That took %.2f seconds' %(time.time()t0) 83 X = domain.vertices 84 HeightQ = domain.quantities['height'].vertex_values 85 VelocityQ = domain.quantities['velocity'].vertex_values 86 x = X.flat 87 z = domain.quantities['elevation'].vertex_values.flat 88 stage=HeightQ.flat+z 77 X = domain.vertices 78 HeightQ = domain.quantities['height'].vertex_values 79 VelocityQ = domain.quantities['velocity'].vertex_values 80 Z = domain.quantities['elevation'].vertex_values 81 Stage = HeightQ+Z 89 82 90 83 from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot … … 94 87 plot1 = subplot(211) 95 88 96 plot( x,z,x,stage)89 plot(X.flat,Z.flat,X.flat,Stage.flat) 97 90 98 91 plot1.set_ylim([1,11]) … … 101 94 legend(('Analytical Solution', 'Numerical Solution'), 102 95 'upper right', shadow=True) 96 103 97 plot2 = subplot(212) 104 plot( x,VelocityQ.flat)98 plot(X.flat,VelocityQ.flat) 105 99 plot2.set_ylim([10,10]) 106 100 
anuga_work/development/2010projects/anuga_1d/channel/test_padarn.py
r7839 r7842 2 2 import random 3 3 from math import sqrt, pow, pi ,sin, cos 4 from channel_domain import * 5 from Numeric import allclose, array, zeros, ones, Float, take, sqrt 6 from config import g, epsilon 4 import numpy 5 6 from anuga_1d.channel.channel_domain import * 7 8 from anuga_1d.config import g, epsilon 7 9 8 10 … … 17 19 omega = sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation 18 20 t=0.0 19 y = zeros(len(x),Float)21 y = numpy.zeros(len(x),numpy.float) 20 22 for i in range(len(x)): 21 23 #y[i] = z_infty+2*A0*z_infty/L_x*cos(omega*t)*(x[i]/L_x0.5*A0/(L_x)*cos(omega*t)) … … 26 28 N = len(x) 27 29 z_infty = 10.0 28 z = zeros(N,Float)30 z = numpy.zeros(N,numpy.float) 29 31 L_x = 2500.0 30 32 A0 = 0.5*L_x … … 35 37 36 38 def initial_area(x): 37 y = zeros(len(x),Float)39 y = numpy.zeros(len(x),numpy.float) 38 40 for i in range(len(x)): 39 41 y[i]=initial_stage([x[i]])bed([x[i]]) … … 61 63 print "Evaluating domain with %d cells" %N 62 64 cell_len = 4*L/N # Origin = 0.0 63 points = zeros(N+1,Float)65 points = numpy.zeros(N+1,numpy.float) 64 66 for i in range(N+1): 65 67 points[i] = 2*L +i*cell_len … … 95 97 96 98 N = float(N) 97 HeightC = domain.quantities['height'].centroid_values99 HeightC = domain.quantities['height'].centroid_values 98 100 DischargeC = domain.quantities['discharge'].centroid_values 99 C = domain.centroids101 C = domain.centroids 100 102 print 'That took %.2f seconds' %(time.time()t0) 101 X = domain.vertices 102 HeightQ = domain.quantities['area'].vertex_values 103 VelocityQ = domain.quantities['velocity'].vertex_values 104 x = X.flat 105 z = domain.quantities['elevation'].vertex_values.flat 106 stage=HeightQ.flat+z 103 X = domain.vertices 104 HeightQ = domain.quantities['area'].vertex_values 105 VelocityQ = domain.quantities['velocity'].vertex_values 106 Z = domain.quantities['elevation'].vertex_values 107 Stage = HeightQ + Z 107 108 108 109 … … 113 114 114 115 plot1 = subplot(211) 115 116 plot(x,z,x,stage) 116 plot(X.flat,Z.flat, X.flat,Stage.flat) 117 117 118 118 plot1.set_ylim([1,35]) … … 122 122 'upper right', shadow=True) 123 123 plot2 = subplot(212) 124 plot( x,VelocityQ.flat)124 plot(X.flat,VelocityQ.flat) 125 125 plot2.set_ylim([10,10]) 126 126 
anuga_work/development/2010projects/anuga_1d/generic/generic_domain.py
r7839 r7842 104 104 self.build_surrogate_neighbour_structure() 105 105 106 #Build boundary dictionary mapping (id, edge) to symbolic tags107 106 #Build boundary dictionary mapping (id, vertex) to symbolic tags 108 107 self.build_boundary_dictionary(boundary) … … 745 744 msg += 'in the supplied dictionary.\n' 746 745 msg += 'The tags are: %s' %self.get_boundary_tags() 747 raise msg746 raise Exception, msg 748 747 749 748 … … 1247 1246 #print 'name %s j = %f \n'%(name,j) 1248 1247 Q = self.quantities[name] 1248 1249 1249 Q.boundary_values[i] = q[j] 1250 1250 #print 'Q=',Q
Note: See TracChangeset
for help on using the changeset viewer.