- Timestamp:
- Sep 5, 2007, 4:42:58 PM (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py
r4702 r4704 26 26 from anuga.abstract_2d_finite_volumes.region\ 27 27 import Set_region as region_set_region 28 29 from anuga.utilities.polygon import inside_polygon 30 from anuga.abstract_2d_finite_volumes.util import get_textual_float 28 31 29 32 import types … … 205 208 self.monitor_polygon = None 206 209 self.monitor_time_interval = None 207 210 self.monitor_indices = None 208 211 209 212 # Checkpointing and storage … … 579 582 self.quantities_to_be_monitored = None 580 583 self.monitor_polygon = None 581 self.monitor_time_interval = None 584 self.monitor_time_interval = None 585 self.monitor_indices = None 582 586 return 583 587 … … 596 600 apply_expression_to_dictionary(quantity_name, self.quantities) 597 601 598 # Initialise extrema 599 self.quantities_to_be_monitored[quantity_name] = [None, None] 602 # Initialise extrema information 603 info_block = {'min': None, # Min value 604 'max': None, # Max value 605 'min_location': None, # Argmin (x, y) 606 'max_location': None, # Argmax (x, y) 607 'min_time': None, # Argmin (t) 608 'max_time': None} # Argmax (t) 600 609 610 self.quantities_to_be_monitored[quantity_name] = info_block 611 612 601 613 602 614 if polygon is not None: 603 # FIXME Check input 604 pass 615 # Check input 616 if isinstance(polygon, basestring): 617 618 # Check if multiple quantities were accidentally 619 # given as separate argument rather than a list. 620 msg = 'Multiple quantities must be specified in a list. ' 621 msg += 'Not as multiple arguments. ' 622 msg += 'I got "%s" as a second argument' %polygon 623 624 if polygon in self.quantities: 625 raise Exception, msg 626 627 try: 628 apply_expression_to_dictionary(polygon, self.quantities) 629 except: 630 # At least polygon wasn't an expression involving quantitites 631 pass 632 else: 633 raise Exception, msg 634 635 # In any case, we don't allow polygon to be a string 636 msg = 'argument "polygon" must not be a string: ' 637 msg += 'I got polygon=\'%s\' ' %polygon 638 raise Exception, msg 639 640 641 # Get indices for centroids that are inside polygon 642 points = self.get_centroid_coordinates(absolute=True) 643 self.monitor_indices = inside_polygon(points, polygon) 644 605 645 606 646 if time_interval is not None: 607 # FIXME Check input 608 pass 609 647 assert len(time_interval) == 2 610 648 611 649 … … 758 796 """ 759 797 760 # Input checks798 # Input checks 761 799 import types, string 762 800 … … 780 818 assert type(tags) == types.ListType, msg 781 819 782 # Determine width of longest quantity name (for cosmetic purposes)820 # Determine width of longest quantity name (for cosmetic purposes) 783 821 maxwidth = 0 784 822 for name in quantities: … … 787 825 maxwidth = w 788 826 789 # Output stats827 # Output stats 790 828 msg = 'Boundary values at time %.4f:\n' %self.time 791 829 for tag in tags: … … 795 833 q = self.quantities[name] 796 834 797 # Find range of boundary values for tag and q835 # Find range of boundary values for tag and q 798 836 maxval = minval = None 799 837 for i, ((vol_id, edge_id), B) in\ … … 815 853 816 854 817 818 def quantity_statistics(self): 855 def update_extrema(self): 856 """Update extrema if requested by set_quantities_to_be_monitored. 857 This data is used for reporting e.g. by running 858 print domain.quantity_statistics() 859 and may also stored in output files (see data_manager in shallow_water) 860 """ 861 862 if self.quantities_to_be_monitored is None: 863 return 864 865 # Observe time interval restriction if any 866 if self.monitor_time_interval is not None and\ 867 (self.time < self.monitor_time_interval[0] or\ 868 self.time > self.monitor_time_interval[1]): 869 return 870 871 872 # Update extrema for each specified quantity subject to 873 # polygon restriction (via monitor_indices). 874 for quantity_name in self.quantities_to_be_monitored: 875 876 if quantity_name in self.quantities: 877 Q = self.get_quantity(quantity_name) 878 else: 879 Q = self.create_quantity_from_expression(quantity_name) 880 881 info_block = self.quantities_to_be_monitored[quantity_name] 882 883 # Update maximum (n > None is always True) 884 maxval = Q.get_maximum_value(self.monitor_indices) 885 if maxval > info_block['max']: 886 info_block['max'] = maxval 887 maxloc = Q.get_maximum_location() 888 info_block['max_location'] = maxloc 889 info_block['max_time'] = self.time 890 891 892 # Update minimum (n < None is always False) 893 minval = Q.get_minimum_value(self.monitor_indices) 894 if info_block['min'] is None or\ 895 minval < info_block['min']: 896 info_block['min'] = minval 897 minloc = Q.get_minimum_location() 898 info_block['min_location'] = minloc 899 info_block['min_time'] = self.time 900 901 902 903 def quantity_statistics(self, precision = '%.4f'): 819 904 """Return string with statistics about quantities for printing or logging 820 905 … … 825 910 """ 826 911 827 pass 828 829 912 maxlen = 128 # Max length of polygon string representation 913 914 # Output stats 915 msg = 'Monitored quantities at time %.4f:\n' %self.time 916 if self.monitor_polygon is not None: 917 p_str = str(self.monitor_polygon) 918 msg += '- Restricted by polygon: %s' %p_str[:maxlen] 919 if len(p_str) >= maxlen: 920 msg += '...\n' 921 else: 922 msg += '\n' 923 924 925 if self.monitor_time_interval is not None: 926 msg += '- Restricted by time interval: %s\n' %str(self.monitor_time_interval) 927 time_interval_start = self.monitor_time_interval[0] 928 else: 929 time_interval_start = 0.0 930 931 932 for quantity_name, info in self.quantities_to_be_monitored.items(): 933 msg += ' %s:\n' %quantity_name 934 935 msg += ' values since time = %.2f in [%s, %s]\n'\ 936 %(time_interval_start, 937 get_textual_float(info['min'], precision), 938 get_textual_float(info['max'], precision)) 939 940 msg += ' minimum attained at time = %s, location = %s\n'\ 941 %(get_textual_float(info['min_time'], precision), 942 get_textual_float(info['min_location'], precision)) 943 944 945 msg += ' maximum attained at time = %s, location = %s\n'\ 946 %(get_textual_float(info['max_time'], precision), 947 get_textual_float(info['max_location'], precision)) 948 949 950 return msg 830 951 831 952 … … 902 1023 from anuga.config import min_timestep, max_timestep, epsilon 903 1024 904 # FIXME: Maybe lump into a larger check prior to evolving1025 # FIXME: Maybe lump into a larger check prior to evolving 905 1026 msg = 'Boundary tags must be bound to boundary objects before ' 906 1027 msg += 'evolving system, ' … … 920 1041 921 1042 if finaltime is not None and duration is not None: 922 # print 'F', finaltime, duration1043 # print 'F', finaltime, duration 923 1044 msg = 'Only one of finaltime and duration may be specified' 924 1045 raise msg … … 932 1053 933 1054 934 self.yieldtime = 0.0 # Time between 'yields'935 936 # Initialise interval of timestep sizes (for reporting only)1055 self.yieldtime = 0.0 # Track time between 'yields' 1056 1057 # Initialise interval of timestep sizes (for reporting only) 937 1058 self.min_timestep = max_timestep 938 1059 self.max_timestep = min_timestep … … 940 1061 self.number_of_first_order_steps = 0 941 1062 942 # update ghosts1063 # Update ghosts 943 1064 self.update_ghosts() 944 1065 945 # Initial update of vertex and edge values1066 # Initial update of vertex and edge values 946 1067 self.distribute_to_vertices_and_edges() 947 1068 948 #Initial update boundary values 1069 # Update extrema if necessary (for reporting) 1070 self.update_extrema() 1071 1072 # Initial update boundary values 949 1073 self.update_boundary() 950 1074 951 # Or maybe restore from latest checkpoint1075 # Or maybe restore from latest checkpoint 952 1076 if self.checkpoint is True: 953 1077 self.goto_latest_checkpoint() 954 1078 955 1079 if skip_initial_step is False: 956 yield(self.time) # Yield initial values1080 yield(self.time) # Yield initial values 957 1081 958 1082 while True: 959 # Compute fluxes across each element edge1083 # Compute fluxes across each element edge 960 1084 self.compute_fluxes() 961 1085 962 # Update timestep to fit yieldstep and finaltime1086 # Update timestep to fit yieldstep and finaltime 963 1087 self.update_timestep(yieldstep, finaltime) 964 1088 965 # Update conserved quantities1089 # Update conserved quantities 966 1090 self.update_conserved_quantities() 967 1091 968 # update ghosts1092 # Update ghosts 969 1093 self.update_ghosts() 970 1094 971 # Update vertex and edge values1095 # Update vertex and edge values 972 1096 self.distribute_to_vertices_and_edges() 973 1097 974 #Update boundary values 1098 # Update extrema if necessary (for reporting) 1099 self.update_extrema() 1100 1101 # Update boundary values 975 1102 self.update_boundary() 976 1103 977 # Update time1104 # Update time 978 1105 self.time += self.timestep 979 1106 self.yieldtime += self.timestep … … 982 1109 self.number_of_first_order_steps += 1 983 1110 984 # Yield results1111 # Yield results 985 1112 if finaltime is not None and self.time >= finaltime-epsilon: 986 1113 987 1114 if self.time > finaltime: 988 # FIXME (Ole, 30 April 2006): Do we need this check?1115 # FIXME (Ole, 30 April 2006): Do we need this check? 989 1116 print 'WARNING (domain.py): time overshot finaltime. Contact Ole.Nielsen@ga.gov.au' 990 1117 self.time = finaltime
Note: See TracChangeset
for help on using the changeset viewer.