Changeset 242
- Timestamp:
- Aug 30, 2004, 4:55:48 PM (21 years ago)
- Location:
- inundation/ga/storm_surge/pyvolution
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/ga/storm_surge/pyvolution/config.py
r240 r242 62 62 63 63 use_psyco = True #Use psyco optimisations 64 #use_psyco = False #Do not use psyco optimisations64 use_psyco = False #Do not use psyco optimisations 65 65 66 66 -
inundation/ga/storm_surge/pyvolution/domain.py
r240 r242 19 19 20 20 from Numeric import zeros, Float 21 from quantity import Quantity 21 from quantity import Quantity, Conserved_quantity 22 22 23 23 #List of quantity names entering … … 34 34 #Build dictionary of Quantity instances keyed by quantity names 35 35 self.quantities = {} 36 for name in self.conserved_quantities + other_quantities: 36 for name in self.conserved_quantities: 37 self.quantities[name] = Conserved_quantity(self) 38 for name in other_quantities: 37 39 self.quantities[name] = Quantity(self) 38 40 -
inundation/ga/storm_surge/pyvolution/quantity.py
r234 r242 15 15 """ 16 16 17 #FIXME: Make Conserved_quantity a subclass of Quantity18 19 from mesh import Mesh20 17 21 18 class Quantity: … … 23 20 def __init__(self, domain, vertex_values=None): 24 21 22 from mesh import Mesh 25 23 from Numeric import array, zeros, Float 26 24 … … 53 51 self.edge_values = zeros((N, 3), Float) 54 52 55 #Allocate space for boundary values56 L = len(domain.boundary)57 self.boundary_values = zeros(L, Float)58 59 #Allocate space for updates of conserved quantities by60 #flux calculations and forcing functions61 self.explicit_update = zeros(N, Float )62 self.semi_implicit_update = zeros(N, Float )63 64 53 #Intialise centroid and edge_values 65 54 self.interpolate() … … 71 60 """ 72 61 73 #FIXME: Write using vector operations (or in C)74 62 N = self.vertex_values.shape[0] 75 63 for i in range(N): … … 79 67 80 68 self.centroid_values[i] = (v0 + v1 + v2)/3 81 82 #FIXME: This is a duplicate from 'interpolate_from_vertices_to_edges' 83 self.edge_values[i, 0] = 0.5*(v1 + v2) 84 self.edge_values[i, 1] = 0.5*(v0 + v2) 85 self.edge_values[i, 2] = 0.5*(v0 + v1) 86 69 70 self.interpolate_from_vertices_to_edges() 71 72 73 def interpolate_from_vertices_to_edges(self): 74 for k in range(self.vertex_values.shape[0]): 75 q0 = self.vertex_values[k, 0] 76 q1 = self.vertex_values[k, 1] 77 q2 = self.vertex_values[k, 2] 78 79 self.edge_values[k, 0] = 0.5*(q1+q2) 80 self.edge_values[k, 1] = 0.5*(q0+q2) 81 self.edge_values[k, 2] = 0.5*(q0+q1) 82 83 84 85 def set_values(self, X, location='vertices'): 86 """Set values for quantity 87 88 X: Compatible list, Numeric array (see below), constant or function 89 location: Where values are to be stored. 90 Permissible options are: vertices, edges, centroid 91 Default is "vertices" 92 93 In case of location == 'centroid' the dimension values must 94 be a list of a Numerical array of length N, N being the number 95 of elements in the mesh. Otherwise it must be of dimension Nx3 96 97 The values will be stored in elements following their 98 internal ordering. 99 100 If values are described a function, it will be evaluated at specified points 101 102 If selected location is vertices, values for centroid and edges 103 will be assigned interpolated values. 104 In any other case, only values for the specified locations 105 will be assigned and the others will be left undefined. 106 """ 107 108 if location not in ['vertices', 'centroids', 'edges']: 109 msg = 'Invalid location: %s' %location 110 raise msg 111 112 if X is None: 113 msg = 'Given values are None' 114 raise msg 115 116 import types 117 118 if callable(X): 119 #Use function specific method 120 self.set_function_values(X, location) 121 elif type(X) in [types.FloatType, types.IntType, types.LongType]: 122 if location == 'centroids': 123 self.centroid_values[:] = X 124 elif location == 'edges': 125 self.edge_values[:] = X 126 else: 127 self.vertex_values[:] = X 128 129 else: 130 #Use array specific method 131 self.set_array_values(X, location) 132 133 if location == 'vertices': 134 #Intialise centroid and edge_values 135 self.interpolate() 136 137 138 139 def set_function_values(self, f, location='vertices'): 140 """Set values for quantity using specified function 141 142 f: x, y -> z Function where x, y and z are arrays 143 location: Where values are to be stored. 144 Permissible options are: vertices, edges, centroid 145 Default is "vertices" 146 """ 147 148 if location == 'centroids': 149 P = self.domain.centroids 150 self.set_values(f(P[:,0], P[:,1]), location) 151 elif location == 'edges': 152 raise 'Not implemented: %s' %location 153 else: 154 #Vertices 155 P = self.domain.get_vertex_coordinates() 156 for i in range(3): 157 self.vertex_values[:,i] = f(P[:,2*i], P[:,2*i+1]) 158 159 160 def set_array_values(self, values, location='vertices'): 161 """Set values for quantity 162 163 values: Numeric array 164 location: Where values are to be stored. 165 Permissible options are: vertices, edges, centroid 166 Default is "vertices" 167 168 In case of location == 'centroid' the dimension values must 169 be a list of a Numerical array of length N, N being the number 170 of elements in the mesh. Otherwise it must be of dimension Nx3 171 172 The values will be stored in elements following their 173 internal ordering. 174 175 If selected location is vertices, values for centroid and edges 176 will be assigned interpolated values. 177 In any other case, only values for the specified locations 178 will be assigned and the others will be left undefined. 179 """ 180 181 from Numeric import array, Float 182 183 values = array(values).astype(Float) 184 185 N = self.centroid_values.shape[0] 186 187 msg = 'Number of values must match number of elements' 188 assert values.shape[0] == N, msg 189 190 if location == 'centroids': 191 assert len(values.shape) == 1, 'Values array must be 1d' 192 self.centroid_values = values 193 elif location == 'edges': 194 assert len(values.shape) == 2, 'Values array must be 2d' 195 msg = 'Array must be N x 3' 196 assert values.shape[1] == 3, msg 197 198 self.edge_values = values 199 else: 200 assert len(values.shape) == 2, 'Values array must be 2d' 201 msg = 'Array must be N x 3' 202 assert values.shape[1] == 3, msg 203 204 self.vertex_values = values 205 206 207 208 class Conserved_quantity(Quantity): 209 """Class conserved quantity adds to Quantity: 210 211 boundary values, storage and method for updating, and 212 methods for extrapolation from centropid to vertices inluding 213 gradients and limiters 214 """ 215 216 def __init__(self, domain, vertex_values=None): 217 Quantity.__init__(self, domain, vertex_values) 218 219 from Numeric import zeros, Float 220 221 #Allocate space for boundary values 222 L = len(domain.boundary) 223 self.boundary_values = zeros(L, Float) 224 225 #Allocate space for updates of conserved quantities by 226 #flux calculations and forcing functions 227 228 N = domain.number_of_elements 229 self.explicit_update = zeros(N, Float ) 230 self.semi_implicit_update = zeros(N, Float ) 231 87 232 88 233 def update(self, timestep): … … 134 279 #We have zero neighbouring volumes - 135 280 #Fall back to first order scheme 136 137 281 pass 138 282 … … 162 306 a[k] = (y1*q0 - y0*q1)/det 163 307 b[k] = (x0*q1 - x1*q0)/det 164 165 308 166 309 else: … … 288 431 289 432 290 291 def interpolate_from_vertices_to_edges(self):292 #FIXME: Write using vector operations (or in C)293 for k in range(self.vertex_values.shape[0]):294 q0 = self.vertex_values[k, 0]295 q1 = self.vertex_values[k, 1]296 q2 = self.vertex_values[k, 2]297 298 self.edge_values[k, 0] = 0.5*(q1+q2)299 self.edge_values[k, 1] = 0.5*(q0+q2)300 self.edge_values[k, 2] = 0.5*(q0+q1)301 302 303 304 def set_values(self, X, location='vertices'):305 """Set values for quantity306 307 X: Compatible list, Numeric array (see below), constant or function308 location: Where values are to be stored.309 Permissible options are: vertices, edges, centroid310 Default is "vertices"311 312 In case of location == 'centroid' the dimension values must313 be a list of a Numerical array of length N, N being the number314 of elements in the mesh. Otherwise it must be of dimension Nx3315 316 The values will be stored in elements following their317 internal ordering.318 319 If values are described a function, it will be evaluated at specified points320 321 If selected location is vertices, values for centroid and edges322 will be assigned interpolated values.323 In any other case, only values for the specified locations324 will be assigned and the others will be left undefined.325 """326 327 328 #FIXME: Should do location error check here only329 330 if location not in ['vertices', 'centroids', 'edges']:331 msg = 'Invalid location: %s' %location332 raise msg333 334 if X is None:335 msg = 'Given values are None'336 raise msg337 338 import types339 340 if callable(X):341 #Use function specific method342 self.set_function_values(X, location)343 elif type(X) in [types.FloatType, types.IntType, types.LongType]:344 if location == 'centroids':345 self.centroid_values[:] = X346 elif location == 'edges':347 self.edge_values[:] = X348 else:349 self.vertex_values[:] = X350 351 else:352 #Use array specific method353 self.set_array_values(X, location)354 355 if location == 'vertices':356 #Intialise centroid and edge_values357 self.interpolate()358 359 360 361 def set_function_values(self, f, location='vertices'):362 """Set values for quantity using specified function363 364 f: x, y -> z Function where x, y and z are arrays365 location: Where values are to be stored.366 Permissible options are: vertices, edges, centroid367 Default is "vertices"368 """369 370 if location == 'centroids':371 P = self.domain.centroids372 self.set_values(f(P[:,0], P[:,1]), location)373 elif location == 'edges':374 raise 'Not yet implemented: %s' %location375 else:376 #Vertices377 P = self.domain.get_vertex_coordinates()378 for i in range(3):379 self.vertex_values[:,i] = f(P[:,2*i], P[:,2*i+1])380 381 382 def set_array_values(self, values, location='vertices'):383 """Set values for quantity384 385 values: Numeric array386 location: Where values are to be stored.387 Permissible options are: vertices, edges, centroid388 Default is "vertices"389 390 In case of location == 'centroid' the dimension values must391 be a list of a Numerical array of length N, N being the number392 of elements in the mesh. Otherwise it must be of dimension Nx3393 394 The values will be stored in elements following their395 internal ordering.396 397 If selected location is vertices, values for centroid and edges398 will be assigned interpolated values.399 In any other case, only values for the specified locations400 will be assigned and the others will be left undefined.401 """402 403 from Numeric import array, Float404 405 values = array(values).astype(Float)406 407 N = self.centroid_values.shape[0]408 409 msg = 'Number of values must match number of elements'410 assert values.shape[0] == N, msg411 412 if location == 'centroids':413 assert len(values.shape) == 1, 'Values array must be 1d'414 self.centroid_values = values415 elif location == 'edges':416 assert len(values.shape) == 2, 'Values array must be 2d'417 msg = 'Array must be N x 3'418 assert values.shape[1] == 3, msg419 420 self.edge_values = values421 else:422 assert len(values.shape) == 2, 'Values array must be 2d'423 msg = 'Array must be N x 3'424 assert values.shape[1] == 3, msg425 426 self.vertex_values = values427 -
inundation/ga/storm_surge/pyvolution/shallow_water.py
r240 r242 352 352 Q.interpolate_from_vertices_to_edges() 353 353 354 354 355 def protect_against_infinitesimal_heights(domain): 355 356 """Protect against infinitesimal heights and high speeds … … 397 398 #Update 398 399 for k in range(domain.number_of_elements): 399 400 401 400 if domain.newstyle: 402 401 if hc[k] < domain.minimum_allowed_height: -
inundation/ga/storm_surge/pyvolution/test_quantity.py
r212 r242 87 87 88 88 def test_boundary_allocation(self): 89 quantity = Quantity(self.mesh4,89 quantity = Conserved_quantity(self.mesh4, 90 90 [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) 91 91 … … 189 189 190 190 def test_gradient(self): 191 quantity = Quantity(self.mesh4)191 quantity = Conserved_quantity(self.mesh4) 192 192 193 193 #Set up for a gradient of (3,0) at mid triangle … … 218 218 219 219 def test_second_order_extrapolation2(self): 220 quantity = Quantity(self.mesh4)220 quantity = Conserved_quantity(self.mesh4) 221 221 222 222 #Set up for a gradient of (3,1), f(x) = 3x+y … … 327 327 328 328 def test_first_order_extrapolator(self): 329 quantity = Quantity(self.mesh4)329 quantity = Conserved_quantity(self.mesh4) 330 330 331 331 #Test centroids … … 342 342 343 343 def test_second_order_extrapolator(self): 344 quantity = Quantity(self.mesh4)344 quantity = Conserved_quantity(self.mesh4) 345 345 346 346 #Set up for a gradient of (3,0) at mid triangle … … 375 375 376 376 377 def test_limiter(self): 378 quantity = Quantity(self.mesh4)377 def test_limiter(self): 378 quantity = Conserved_quantity(self.mesh4) 379 379 380 380 #Create a deliberate overshoot (e.g. from gradient computation) … … 406 406 407 407 def test_distribute_first_order(self): 408 quantity = Quantity(self.mesh4)408 quantity = Conserved_quantity(self.mesh4) 409 409 410 410 #Test centroids … … 427 427 428 428 def test_update_explicit(self): 429 quantity = Quantity(self.mesh4)429 quantity = Conserved_quantity(self.mesh4) 430 430 431 431 #Test centroids … … 443 443 444 444 def test_update_semi_implicit(self): 445 quantity = Quantity(self.mesh4)445 quantity = Conserved_quantity(self.mesh4) 446 446 447 447 #Test centroids … … 459 459 460 460 def test_both_updates(self): 461 quantity = Quantity(self.mesh4)461 quantity = Conserved_quantity(self.mesh4) 462 462 463 463 #Test centroids
Note: See TracChangeset
for help on using the changeset viewer.