Changeset 4771
- Timestamp:
- Oct 30, 2007, 6:25:23 PM (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py
r4735 r4771 190 190 # Or maybe get rid of order altogether and use beta_w and beta_h 191 191 self.set_default_order(1) 192 #self.default_order = 1193 #self._order_ = self.default_order194 192 195 193 self.smallsteps = 0 … … 207 205 self.finaltime = None 208 206 self.min_timestep = self.max_timestep = 0.0 209 self.starttime = 0 #Physical starttime if any (0 is 1 Jan 1970 00:00:00) 207 self.starttime = 0 # Physical starttime if any 208 # (0 is 1 Jan 1970 00:00:00) 210 209 self.timestep = 0.0 211 210 self.flux_timestep = 0.0 … … 223 222 self.checkpoint = False 224 223 225 # MH310505 To avoid calculating the flux across each edge twice, keep an integer (boolean) array, 226 # to be used during the flux calculation 227 N = len(self) #number_of_triangles 224 # To avoid calculating the flux across each edge twice, 225 # keep an integer (boolean) array, to be used during the flux 226 # calculation 227 N = len(self) # Number_of_triangles 228 228 self.already_computed_flux = zeros((N, 3), Int) 229 229 230 # Storage for maximal speeds computed for each triangle by compute_fluxes 230 # Storage for maximal speeds computed for each triangle by 231 # compute_fluxes 231 232 # This is used for diagnostics only 232 233 self.max_speed = zeros(N, Float) … … 246 247 247 248 248 249 #Public interface to Domain 249 #--------------------------- 250 # Public interface to Domain 251 #--------------------------- 250 252 def get_conserved_quantities(self, vol_id, vertex=None, edge=None): 251 253 """Get conserved quantities at volume vol_id … … 332 334 """ 333 335 334 #FIXME (Ole): Allow new quantities here 335 #from quantity import Quantity, Conserved_quantity 336 #Create appropriate quantity object 337 ##if name in self.conserved_quantities: 338 ## self.quantities[name] = Conserved_quantity(self) 339 ##else: 340 ## self.quantities[name] = Quantity(self) 341 342 343 #Do the expression stuff 336 # Do the expression stuff 344 337 if kwargs.has_key('expression'): 345 338 expression = kwargs['expression'] … … 349 342 kwargs['quantity'] = Q 350 343 351 # Assign values344 # Assign values 352 345 self.quantities[name].set_values(*args, **kwargs) 353 346 … … 463 456 464 457 465 # FIXME: Try to remove the sorting and fix test_mesh.py458 # FIXME (Ole): Try to remove the sorting and fix test_mesh.py 466 459 x = self.boundary.keys() 467 460 x.sort() 468 461 469 # Loop through edges that lie on the boundary and associate them with470 # callable boundary objects depending on their tags462 # Loop through edges that lie on the boundary and associate them with 463 # callable boundary objects depending on their tags 471 464 self.boundary_objects = [] 472 465 for k, (vol_id, edge_id) in enumerate(x): … … 474 467 475 468 if self.boundary_map.has_key(tag): 476 B = self.boundary_map[tag] # Get callable boundary object469 B = self.boundary_map[tag] # Get callable boundary object 477 470 478 471 if B is not None: … … 509 502 can also be passed into this function. 510 503 """ 511 #print "*args", args 512 #print "**kwargs", kwargs 504 513 505 if len(args) == 1: 514 506 self._set_region(*args, **kwargs) 515 507 else: 516 # Assume it is arguments for the region.set_region function508 # Assume it is arguments for the region.set_region function 517 509 func = region_set_region(*args, **kwargs) 518 510 self._set_region(func) … … 628 620 629 621 630 631 #MISC 622 #-------------------------- 623 # Miscellaneous diagnostics 624 #-------------------------- 632 625 def check_integrity(self): 633 626 Mesh.check_integrity(self) … … 638 631 639 632 ##assert hasattr(self, 'boundary_objects') 633 640 634 641 635 def write_time(self, track_speeds=False): … … 646 640 """Return string with time stepping statistics for printing or logging 647 641 648 Optional boolean keyword track_speeds decides whether to report location of 649 smallest timestep as well as a histogram and percentile report. 642 Optional boolean keyword track_speeds decides whether to report 643 location of smallest timestep as well as a histogram and percentile 644 report. 650 645 """ 651 646 … … 672 667 673 668 674 # Setup 10 bins for speed histogram669 # Setup 10 bins for speed histogram 675 670 bins = create_bins(self.max_speed, 10) 676 671 hist = histogram(self.max_speed, bins) 677 672 678 673 msg += '------------------------------------------------\n' 679 msg += ' Speeds in [%f, %f]\n' %(min(self.max_speed), max(self.max_speed)) 674 msg += ' Speeds in [%f, %f]\n' %(min(self.max_speed), 675 max(self.max_speed)) 680 676 msg += ' Histogram:\n' 681 677 … … 684 680 lo = hi 685 681 if i+1 < len(bins): 686 # Open upper interval682 # Open upper interval 687 683 hi = bins[i+1] 688 684 msg += ' [%f, %f[: %d\n' %(lo, hi, count) 689 685 else: 690 # Closed upper interval686 # Closed upper interval 691 687 hi = max(self.max_speed) 692 688 msg += ' [%f, %f]: %d\n' %(lo, hi, count) … … 702 698 lower = min(speed) 703 699 for i, a in enumerate(speed): 704 if i % (N/10) == 0 and i != 0: #For every 10% of the sorted speeds 700 if i % (N/10) == 0 and i != 0: 701 # For every 10% of the sorted speeds 705 702 msg += ' %d speeds in [%f, %f]\n' %(i-k, lower, a) 706 703 lower = a … … 711 708 712 709 713 714 710 715 711 … … 755 751 756 752 Input: 757 quantities: either None, a string or a list of strings naming the quantities to be reported 758 tags: either None, a string or a list of strings naming the tags to be reported 753 quantities: either None, a string or a list of strings naming the 754 quantities to be reported 755 tags: either None, a string or a list of strings naming the 756 tags to be reported 759 757 760 758 … … 767 765 768 766 769 If quantities are specified only report on those. Otherwise take all conserved quantities. 767 If quantities are specified only report on those. Otherwise take all 768 conserved quantities. 770 769 If tags are specified only report on those, otherwise take all tags. 771 770 … … 801 800 maxwidth = w 802 801 803 # Output stat s802 # Output statistics 804 803 msg = 'Boundary values at time %.4f:\n' %self.time 805 804 for tag in tags: … … 860 859 861 860 # Update maximum 862 # (n > None is always True, but we check explicitly because of the epsilon) 861 # (n > None is always True, but we check explicitly because 862 # of the epsilon) 863 863 maxval = Q.get_maximum_value(self.monitor_indices) 864 864 if info_block['max'] is None or\ … … 882 882 883 883 def quantity_statistics(self, precision = '%.4f'): 884 """Return string with statistics about quantities for printing or logging 884 """Return string with statistics about quantities for 885 printing or logging 885 886 886 887 Quantities reported are specified through method … … 892 893 maxlen = 128 # Max length of polygon string representation 893 894 894 # Output stat s895 # Output statistics 895 896 msg = 'Monitored quantities at time %.4f:\n' %self.time 896 897 if self.monitor_polygon is not None: … … 904 905 905 906 if self.monitor_time_interval is not None: 906 msg += '- Restricted by time interval: %s\n' %str(self.monitor_time_interval) 907 msg += '- Restricted by time interval: %s\n'\ 908 %str(self.monitor_time_interval) 907 909 time_interval_start = self.monitor_time_interval[0] 908 910 else: … … 916 918 %(time_interval_start, 917 919 get_textual_float(info['min'], precision), 918 get_textual_float(info['max'], precision)) 920 get_textual_float(info['max'], precision)) 919 921 920 922 msg += ' minimum attained at time = %s, location = %s\n'\ 921 923 %(get_textual_float(info['min_time'], precision), 922 get_textual_float(info['min_location'], precision)) 923 924 get_textual_float(info['min_location'], precision)) 925 924 926 925 927 msg += ' maximum attained at time = %s, location = %s\n'\ 926 928 %(get_textual_float(info['max_time'], precision), 927 get_textual_float(info['max_location'], precision)) 929 get_textual_float(info['max_location'], precision)) 928 930 929 931 … … 969 971 970 972 971 #def set_defaults(self): 972 # """Set default values for uninitialised quantities. 973 # Should be overridden or specialised by specific modules 974 # """# 975 # 976 # for name in self.conserved_quantities + self.other_quantities: 977 # self.set_quantity(name, 0.0) 978 979 980 ########################### 981 #Main components of evolve 973 #-------------------------- 974 # Main components of evolve 975 #-------------------------- 982 976 983 977 def evolve(self, … … 1103 1097 if self.time > finaltime: 1104 1098 # FIXME (Ole, 30 April 2006): Do we need this check? 1105 # Probably not (Ole, 18 September 2008). Now changed to Exception 1106 msg = 'WARNING (domain.py): time overshot finaltime. Contact Ole.Nielsen@ga.gov.au' 1099 # Probably not (Ole, 18 September 2008). Now changed to 1100 # Exception 1101 msg = 'WARNING (domain.py): time overshot finaltime. ' 1102 msg += 'Contact Ole.Nielsen@ga.gov.au' 1107 1103 raise Exception, msg 1108 #self.time = finaltime1109 1104 1110 1105 … … 1331 1326 1332 1327 def backup_conserved_quantities(self): 1333 N = len(self) # number_of_triangles1334 1335 # backup conserved_quantities centroid values1328 N = len(self) # Number_of_triangles 1329 1330 # Backup conserved_quantities centroid values 1336 1331 for name in self.conserved_quantities: 1337 1332 Q = self.quantities[name] … … 1341 1336 N = len(self) #number_of_triangles 1342 1337 1343 # backup conserved_quantities centroid values1338 # Backup conserved_quantities centroid values 1344 1339 for name in self.conserved_quantities: 1345 1340 Q = self.quantities[name] … … 1356 1351 """ 1357 1352 1358 # FIXME: Update only those that change (if that can be worked out)1359 # FIXME: Boundary objects should not include ghost nodes.1353 # FIXME: Update only those that change (if that can be worked out) 1354 # FIXME: Boundary objects should not include ghost nodes. 1360 1355 for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects): 1361 1356 if B is None: … … 1403 1398 for i in range(self.number_of_full_triangles): 1404 1399 if self.max_speed[i] > bins[-1]: 1405 msg = 'Time=%f: Ignoring isolated high speed triangle ' %self.time 1400 msg = 'Time=%f: Ignoring isolated high ' %self.time 1401 msg += 'speed triangle ' 1406 1402 msg += '#%d of %d with max speed=%f'\ 1407 %(i, self.number_of_full_triangles, self.max_speed[i])1408 #print msg1403 %(i, self.number_of_full_triangles, 1404 self.max_speed[i]) 1409 1405 1410 # print 'Found offending triangle', i, self.max_speed[i] 1406 # print 'Found offending triangle', i, 1407 # self.max_speed[i] 1411 1408 self.get_quantity('xmomentum').set_values(0.0, indices=[i]) 1412 1409 self.get_quantity('ymomentum').set_values(0.0, indices=[i]) … … 1423 1420 timestep = min(self.CFL*self.flux_timestep, max_timestep) 1424 1421 1425 # Record maximal and minimal values of timestep for reporting1422 # Record maximal and minimal values of timestep for reporting 1426 1423 self.max_timestep = max(timestep, self.max_timestep) 1427 1424 self.min_timestep = min(timestep, self.min_timestep) … … 1429 1426 1430 1427 1431 # Protect against degenerate time steps1428 # Protect against degenerate time steps 1432 1429 if timestep < min_timestep: 1433 1430 1434 # Number of consecutive small steps taken b4 taking action1431 # Number of consecutive small steps taken b4 taking action 1435 1432 self.smallsteps += 1 1436 1433 1437 1434 if self.smallsteps > self.max_smallsteps: 1438 self.smallsteps = 0 # Reset1435 self.smallsteps = 0 # Reset 1439 1436 1440 1437 if self._order_ == 1: … … 1444 1441 %self.max_smallsteps 1445 1442 print msg 1446 timestep = min_timestep #Try enforcing min_step 1447 1443 timestep = min_timestep # Try enforcing min_step 1448 1444 1449 1445 #print self.timestepping_statistics(track_speeds=True) 1450 1451 1446 1452 1447 raise Exception, msg 1453 1448 else: 1454 # Try to overcome situation by switching to 1 order1449 # Try to overcome situation by switching to 1 order 1455 1450 self._order_ = 1 1456 1451 … … 1461 1456 1462 1457 1463 # Ensure that final time is not exceeded1458 # Ensure that final time is not exceeded 1464 1459 if finaltime is not None and self.time + timestep > finaltime : 1465 1460 timestep = finaltime-self.time 1466 1461 1467 # Ensure that model time is aligned with yieldsteps1462 # Ensure that model time is aligned with yieldsteps 1468 1463 if self.yieldtime + timestep > yieldstep: 1469 1464 timestep = yieldstep-self.yieldtime … … 1491 1486 from Numeric import ones, sum, equal, Float 1492 1487 1493 N = len(self) # number_of_triangles1488 N = len(self) # Number_of_triangles 1494 1489 d = len(self.conserved_quantities) 1495 1490 1496 1491 timestep = self.timestep 1497 1492 1498 # Compute forcing terms1493 # Compute forcing terms 1499 1494 self.compute_forcing_terms() 1500 1495 1501 # Update conserved_quantities1496 # Update conserved_quantities 1502 1497 for name in self.conserved_quantities: 1503 1498 Q = self.quantities[name] 1504 1499 Q.update(timestep) 1505 1500 1506 #Clean up 1507 #Note that Q.explicit_update is reset by compute_fluxes 1508 1509 #MH090605 commented out the following since semi_implicit_update is now re-initialized 1510 #at the end of the _update function in quantity_ext.c (This is called by the 1511 #preceeding Q.update(timestep) statement above). 1512 #For run_profile.py with N=128, the time of update_conserved_quantities is cut from 14.00 secs 1513 #to 8.35 secs 1514 1515 #Q.semi_implicit_update[:] = 0.0 1501 # Note that Q.explicit_update is reset by compute_fluxes 1502 # Where is Q.semi_implicit_update reset? 1503 1516 1504 1517 1505 def update_ghosts(self): … … 1551 1539 1552 1540 1553 ############################################## 1554 #Initialise module 1555 1556 #Optimisation with psyco 1541 #------------------ 1542 # Initialise module 1543 #------------------ 1544 1545 # Optimisation with psyco 1557 1546 from anuga.config import use_psyco 1558 1547 if use_psyco: … … 1563 1552 if os.name == 'posix' and os.uname()[4] == 'x86_64': 1564 1553 pass 1565 # Psyco isn't supported on 64 bit systems, but it doesn't matter1554 # Psyco isn't supported on 64 bit systems, but it doesn't matter 1566 1555 else: 1567 1556 msg = 'WARNING: psyco (speedup) could not import'+\ … … 1570 1559 else: 1571 1560 psyco.bind(Domain.update_boundary) 1572 #psyco.bind(Domain.update_timestep) #Not worth it1561 #psyco.bind(Domain.update_timestep) # Not worth it 1573 1562 psyco.bind(Domain.update_conserved_quantities) 1574 1563 psyco.bind(Domain.distribute_to_vertices_and_edges)
Note: See TracChangeset
for help on using the changeset viewer.
