Changeset 6051
- Timestamp:
- Dec 9, 2008, 6:44:58 PM (16 years ago)
- Location:
- anuga_core/source/anuga
- Files:
-
- 4 edited
-
abstract_2d_finite_volumes/domain.py (modified) (5 diffs)
-
culvert_flows/Test_Culvert_Flat_Water_Lev.py (modified) (3 diffs)
-
culvert_flows/culvert_class.py (modified) (5 diffs)
-
shallow_water/test_shallow_water_domain.py (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py
r5961 r6051 122 122 123 123 124 # Build dictionary of Quantity instances keyed by quantity names124 # Build dictionary of Quantity instances keyed by quantity names 125 125 self.quantities = {} 126 126 127 # FIXME: remove later - maybe OK, though....127 # FIXME: remove later - maybe OK, though.... 128 128 for name in self.conserved_quantities: 129 129 self.quantities[name] = Quantity(self) … … 131 131 self.quantities[name] = Quantity(self) 132 132 133 # Create an empty list for explicit forcing terms133 # Create an empty list for explicit forcing terms 134 134 self.forcing_terms = [] 135 135 136 # Setup the ghost cell communication136 # Setup the ghost cell communication 137 137 if full_send_dict is None: 138 138 self.full_send_dict = {} … … 1219 1219 # Update ghosts 1220 1220 self.update_ghosts() 1221 1222 1221 1223 1222 # Update time … … 1560 1559 1561 1560 1562 1563 1561 def update_conserved_quantities(self): 1564 1562 """Update vectors of conserved quantities using previously … … 1582 1580 1583 1581 # Note that Q.explicit_update is reset by compute_fluxes 1584 # Where is Q.semi_implicit_update reset? 1582 # Where is Q.semi_implicit_update reset? 1583 # It is reset in quantity_ext.c 1585 1584 1586 1585 -
anuga_core/source/anuga/culvert_flows/Test_Culvert_Flat_Water_Lev.py
r5868 r6051 138 138 culvert_routine=boyd_generalised_culvert_model, 139 139 number_of_barrels=1, 140 update_interval=2, 140 141 verbose=True) 141 142 … … 158 159 #------------------------------------------------------------------------------ 159 160 160 for t in domain.evolve(yieldstep = 0.01, finaltime = 45): 161 if int(domain.time*100) % 100 == 0: 162 domain.write_time() 161 #for t in domain.evolve(yieldstep = 1, finaltime = 25): 162 # print domain.timestepping_statistics() 163 163 164 164 … … 170 170 t0 = time.time() 171 171 172 s = 'for t in domain.evolve(yieldstep = 0.5, finaltime = 2): domain.write_time()'172 s = 'for t in domain.evolve(yieldstep = 1, finaltime = 25): domain.write_time()' 173 173 174 174 import profile, pstats -
anuga_core/source/anuga/culvert_flows/culvert_class.py
r5777 r6051 64 64 culvert_routine=None, 65 65 number_of_barrels=1, 66 update_interval=None, 66 67 verbose=False): 67 68 … … 191 192 self.verbose = verbose 192 193 self.last_time = 0.0 194 self.last_update = 0.0 # For use with update_interval 195 self.update_interval = update_interval 193 196 194 197 … … 232 235 log_to_file(self.log_filename, s) 233 236 237 234 238 def __call__(self, domain): 235 239 from anuga.utilities.numerical_tools import mean … … 244 248 # Time stuff 245 249 time = domain.get_time() 246 delta_t = time-self.last_time 247 s = '\nTime = %.2f, delta_t = %f' %(time, delta_t) 248 log_to_file(log_filename, s) 249 250 msg = 'Time did not advance' 251 if time > 0.0: assert delta_t > 0.0, msg 252 253 254 # Get average water depths at each opening 255 openings = self.openings # There are two Opening [0] and [1] 256 for i, opening in enumerate(openings): 257 dq = domain.quantities 250 251 252 update = False 253 if self.update_interval is None: 254 update = True 255 else: 256 if time - self.last_update > self.update_interval or time == 0.0: 257 update = True 258 259 #print 'call', time, time - self.last_update 260 261 262 if update is True: 263 #print 'Updating', time, time - self.last_update 264 self.last_update = time 265 266 delta_t = time-self.last_time 267 s = '\nTime = %.2f, delta_t = %f' %(time, delta_t) 268 log_to_file(log_filename, s) 269 270 msg = 'Time did not advance' 271 if time > 0.0: assert delta_t > 0.0, msg 272 273 274 # Get average water depths at each opening 275 openings = self.openings # There are two Opening [0] and [1] 276 for i, opening in enumerate(openings): 277 dq = domain.quantities 278 279 # Compute mean values of selected quantitites in the 280 # exchange area in front of the culvert 281 # Stage and velocity comes from enquiry area 282 # and elevation from exchange area 283 284 stage = dq['stage'].get_values(location='centroids', 285 indices=opening.exchange_indices) 286 w = mean(stage) # Average stage 287 288 # Use invert level instead of elevation if specified 289 invert_level = self.invert_levels[i] 290 if invert_level is not None: 291 z = invert_level 292 else: 293 elevation = dq['elevation'].get_values(location='centroids', 294 indices=opening.exchange_indices) 295 z = mean(elevation) # Average elevation 296 297 # Estimated depth above the culvert inlet 298 d = w - z # Used for calculations involving depth 299 if d < 0.0: 300 # This is possible since w and z are taken at different locations 301 #msg = 'D < 0.0: %f' %d 302 #raise msg 303 d = 0.0 304 305 306 # Ratio of depth to culvert height. 307 # If ratio > 1 then culvert is running full 308 if self.culvert_type == 'circle': 309 ratio = d/self.diameter 310 else: 311 ratio = d/self.height 312 opening.ratio = ratio 313 314 315 # Average measures of energy in front of this opening 316 polyline = self.enquiry_polylines[i] 317 #print 't = %.4f, opening=%d,' %(domain.time, i), 318 opening.total_energy = domain.get_energy_through_cross_section(polyline, 319 kind='total') 320 #print 'Et = %.3f m' %opening.total_energy 321 322 # Store current average stage and depth with each opening object 323 opening.depth = d 324 opening.depth_trigger = d # Use this for now 325 opening.stage = w 326 opening.elevation = z 327 328 329 ################# End of the FOR loop ################################################ 330 331 # We now need to deal with each opening individually 332 333 # Determine flow direction based on total energy difference 334 delta_Et = openings[0].total_energy - openings[1].total_energy 335 336 if delta_Et > 0: 337 #print 'Flow U/S ---> D/S' 338 inlet = openings[0] 339 outlet = openings[1] 340 341 inlet.momentum = self.opening_momentum[0] 342 outlet.momentum = self.opening_momentum[1] 343 344 else: 345 #print 'Flow D/S ---> U/S' 346 inlet = openings[1] 347 outlet = openings[0] 348 349 inlet.momentum = self.opening_momentum[1] 350 outlet.momentum = self.opening_momentum[0] 351 352 delta_Et = -delta_Et 353 354 self.inlet = inlet 355 self.outlet = outlet 356 357 msg = 'Total energy difference is negative' 358 assert delta_Et > 0.0, msg 359 360 delta_h = inlet.stage - outlet.stage 361 delta_z = inlet.elevation - outlet.elevation 362 culvert_slope = (delta_z/self.length) 363 364 if culvert_slope < 0.0: 365 # Adverse gradient - flow is running uphill 366 # Flow will be purely controlled by uphill outlet face 367 if self.verbose is True: 368 print 'WARNING: Flow is running uphill. Watch Out!', inlet.elevation, outlet.elevation 369 370 371 s = 'Delta total energy = %.3f' %(delta_Et) 372 log_to_file(log_filename, s) 373 374 375 # Calculate discharge for one barrel and set inlet.rate and outlet.rate 376 Q, barrel_velocity, culvert_outlet_depth = self.culvert_routine(self, inlet, outlet, delta_Et, g) 377 378 # Adjust discharge for multiple barrels 379 Q *= self.number_of_barrels 380 381 # Compute barrel momentum 382 barrel_momentum = barrel_velocity*culvert_outlet_depth 383 384 s = 'Barrel velocity = %f' %barrel_velocity 385 log_to_file(log_filename, s) 386 387 # Compute momentum vector at outlet 388 outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum 389 390 s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y) 391 log_to_file(log_filename, s) 392 393 # Log timeseries to file 394 fid = open(self.timeseries_filename, 'a') 395 fid.write('%f, %f, %f, %f, %f\n'\ 396 %(time, 397 openings[0].total_energy, 398 openings[1].total_energy, 399 barrel_velocity, 400 Q)) 401 fid.close() 402 403 # Update momentum 404 delta_t = time - self.last_time 405 if delta_t > 0.0: 406 xmomentum_rate = outlet_mom_x - outlet.momentum[0].value 407 xmomentum_rate /= delta_t 408 409 ymomentum_rate = outlet_mom_y - outlet.momentum[1].value 410 ymomentum_rate /= delta_t 411 412 s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate) 413 log_to_file(log_filename, s) 414 else: 415 xmomentum_rate = ymomentum_rate = 0.0 416 417 418 # Set momentum rates for outlet jet 419 outlet.momentum[0].rate = xmomentum_rate 420 outlet.momentum[1].rate = ymomentum_rate 421 422 # Remember this value for next step (IMPORTANT) 423 outlet.momentum[0].value = outlet_mom_x 424 outlet.momentum[1].value = outlet_mom_y 425 426 if int(domain.time*100) % 100 == 0: 427 s = 'T=%.5f, Culvert Discharge = %.3f f'\ 428 %(time, Q) 429 s += ' Depth= %0.3f Momentum = (%0.3f, %0.3f)'\ 430 %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y) 431 s += ' Momentum rate: (%.4f, %.4f)'\ 432 %(xmomentum_rate, ymomentum_rate) 433 s+='Outlet Vel= %.3f'\ 434 %(barrel_velocity) 435 log_to_file(log_filename, s) 258 436 259 # Compute mean values of selected quantitites in the exchange area in front of the culvert 260 # Stage and velocity comes from enquiry area and elevation from exchange area 261 262 stage = dq['stage'].get_values(location='centroids', 263 indices=opening.exchange_indices) 264 w = mean(stage) # Average stage 265 266 # Use invert level instead of elevation if specified 267 invert_level = self.invert_levels[i] 268 if invert_level is not None: 269 z = invert_level 270 else: 271 elevation = dq['elevation'].get_values(location='centroids', indices=opening.exchange_indices) 272 z = mean(elevation) # Average elevation 273 274 # Estimated depth above the culvert inlet 275 d = w - z # Used for calculations involving depth 276 if d < 0.0: 277 # This is possible since w and z are taken at different locations 278 #msg = 'D < 0.0: %f' %d 279 #raise msg 280 d = 0.0 281 282 283 # Ratio of depth to culvert height. 284 # If ratio > 1 then culvert is running full 285 if self.culvert_type == 'circle': 286 ratio = d/self.diameter 287 else: 288 ratio = d/self.height 289 opening.ratio = ratio 290 291 292 # Average measures of energy in front of this opening 293 polyline = self.enquiry_polylines[i] 294 #print 't = %.4f, opening=%d,' %(domain.time, i), 295 opening.total_energy = domain.get_energy_through_cross_section(polyline, 296 kind='total') 297 #print 'Et = %.3f m' %opening.total_energy 298 299 # Store current average stage and depth with each opening object 300 opening.depth = d 301 opening.depth_trigger = d # Use this for now 302 opening.stage = w 303 opening.elevation = z 304 305 306 ################# End of the FOR loop ################################################ 307 308 309 # We now need to deal with each opening individually 310 311 # Determine flow direction based on total energy difference 312 delta_Et = openings[0].total_energy - openings[1].total_energy 313 314 if delta_Et > 0: 315 #print 'Flow U/S ---> D/S' 316 inlet = openings[0] 317 outlet = openings[1] 318 319 inlet.momentum = self.opening_momentum[0] 320 outlet.momentum = self.opening_momentum[1] 321 322 else: 323 #print 'Flow D/S ---> U/S' 324 inlet = openings[1] 325 outlet = openings[0] 326 327 inlet.momentum = self.opening_momentum[1] 328 outlet.momentum = self.opening_momentum[0] 329 330 delta_Et = -delta_Et 331 332 msg = 'Total energy difference is negative' 333 assert delta_Et > 0.0, msg 334 335 delta_h = inlet.stage - outlet.stage 336 delta_z = inlet.elevation - outlet.elevation 337 culvert_slope = (delta_z/self.length) 338 339 if culvert_slope < 0.0: 340 # Adverse gradient - flow is running uphill 341 # Flow will be purely controlled by uphill outlet face 342 print 'WARNING: Flow is running uphill. Watch Out!', inlet.elevation, outlet.elevation 343 344 345 s = 'Delta total energy = %.3f' %(delta_Et) 346 log_to_file(log_filename, s) 347 348 349 # Calculate discharge for one barrel 350 Q, barrel_velocity, culvert_outlet_depth = self.culvert_routine(self, inlet, outlet, delta_Et, g) 351 352 # Adjust discharge for multiple barrels 353 Q *= self.number_of_barrels 354 355 # Compute barrel momentum 356 barrel_momentum = barrel_velocity*culvert_outlet_depth 357 358 s = 'Barrel velocity = %f' %barrel_velocity 359 log_to_file(log_filename, s) 360 361 # Compute momentum vector at outlet 362 outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum 363 364 s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y) 365 log_to_file(log_filename, s) 366 367 # Log timeseries to file 368 fid = open(self.timeseries_filename, 'a') 369 fid.write('%f, %f, %f, %f, %f\n'\ 370 %(time, 371 openings[0].total_energy, 372 openings[1].total_energy, 373 barrel_velocity, 374 Q)) 375 fid.close() 376 377 # Update momentum 378 delta_t = time - self.last_time 379 if delta_t > 0.0: 380 xmomentum_rate = outlet_mom_x - outlet.momentum[0].value 381 xmomentum_rate /= delta_t 382 383 ymomentum_rate = outlet_mom_y - outlet.momentum[1].value 384 ymomentum_rate /= delta_t 385 386 s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate) 387 log_to_file(log_filename, s) 388 else: 389 xmomentum_rate = ymomentum_rate = 0.0 390 391 392 # Set momentum rates for outlet jet 393 outlet.momentum[0].rate = xmomentum_rate 394 outlet.momentum[1].rate = ymomentum_rate 395 396 # Remember this value for next step (IMPORTANT) 397 outlet.momentum[0].value = outlet_mom_x 398 outlet.momentum[1].value = outlet_mom_y 399 400 if int(domain.time*100) % 100 == 0: 401 s = 'T=%.5f, Culvert Discharge = %.3f f'\ 402 %(time, Q) 403 s += ' Depth= %0.3f Momentum = (%0.3f, %0.3f)'\ 404 %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y) 405 s += ' Momentum rate: (%.4f, %.4f)'\ 406 %(xmomentum_rate, ymomentum_rate) 407 s+='Outlet Vel= %.3f'\ 408 %(barrel_velocity) 409 log_to_file(log_filename, s) 410 411 412 413 414 437 438 439 415 440 # Execute flow term for each opening 416 441 # This is where Inflow objects are evaluated and update the domain … … 420 445 # Execute momentum terms 421 446 # This is where Inflow objects are evaluated and update the domain 422 outlet.momentum[0](domain)423 outlet.momentum[1](domain)447 self.outlet.momentum[0](domain) 448 self.outlet.momentum[1](domain) 424 449 425 # Store value of time 450 # Store value of time #FIXME(Ole): Maybe only every time we update 426 451 self.last_time = time 427 452 -
anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py
r6045 r6051 2827 2827 2828 2828 2829 2830 2831 2832 def test_rainfall_forcing_with_evolve(self): 2833 """test_rainfall_forcing_with_evolve 2834 2835 Test how forcing terms are called within evolve 2836 """ 2837 2838 # FIXME(Ole): This test is just to experiment 2839 2840 a = [0.0, 0.0] 2841 b = [0.0, 2.0] 2842 c = [2.0, 0.0] 2843 d = [0.0, 4.0] 2844 e = [2.0, 2.0] 2845 f = [4.0, 0.0] 2846 2847 points = [a, b, c, d, e, f] 2848 #bac, bce, ecf, dbe 2849 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2850 2851 2852 domain = Domain(points, vertices) 2853 2854 #Flat surface with 1m of water 2855 domain.set_quantity('elevation', 0) 2856 domain.set_quantity('stage', 1.0) 2857 domain.set_quantity('friction', 0) 2858 2859 Br = Reflective_boundary(domain) 2860 domain.set_boundary({'exterior': Br}) 2861 2862 # Setup only one forcing term, time dependent rainfall that expires at t==20 2863 from anuga.fit_interpolate.interpolate import Modeltime_too_late 2864 def main_rate(t): 2865 if t > 20: 2866 msg = 'Model time exceeded.' 2867 raise Modeltime_too_late, msg 2868 else: 2869 return 3*t + 7 2870 2871 domain.forcing_terms = [] 2872 R = Rainfall(domain, 2873 rate=main_rate, 2874 polygon=[[1,1], [2,1], [2,2], [1,2]], 2875 default_rate=5.0) 2876 2877 assert allclose(R.exchange_area, 1) 2878 2879 domain.forcing_terms.append(R) 2880 2881 for t in domain.evolve(yieldstep=1, finaltime=25): 2882 pass 2883 2884 #print t, domain.quantities['stage'].explicit_update, (3*t+7)/1000 2885 2886 #FIXME(Ole): A test here is hard because explicit_update also 2887 # receives updates from the flux calculation. 2888 2889 2890 2829 2891 2830 2892 def test_inflow_using_circle(self): … … 6385 6447 if __name__ == "__main__": 6386 6448 6387 suite = unittest.makeSuite(Test_Shallow_Water,'test') 6449 suite = unittest.makeSuite(Test_Shallow_Water, 'test') 6450 #suite = unittest.makeSuite(Test_Shallow_Water, 'test_rainfall_forcing_with_evolve') 6388 6451 #suite = unittest.makeSuite(Test_Shallow_Water,'test_get_energy_through_cross_section_with_g') 6389 6452 #suite = unittest.makeSuite(Test_Shallow_Water,'test_fitting_using_shallow_water_domain')
Note: See TracChangeset
for help on using the changeset viewer.
