Changeset 5738
- Timestamp:
- Sep 5, 2008, 3:06:01 PM (17 years ago)
- Location:
- anuga_work/development/anuga_1d
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/development/anuga_1d/compile_all.py
r5728 r5738 6 6 7 7 anuga_dir = '..' + os.sep + '..' + os.sep + '..' 8 #os.chdir(anuga_dir)9 print os.getcwd()10 11 8 utilities_dir = anuga_dir + os.sep + 'anuga_core' + os.sep + 'source' + os.sep + 'anuga' + os.sep + 'utilities' 12 9 13 print anuga_dir14 print utilities_dir15 10 execfile( utilities_dir + os.sep + 'compile.py') 16 11 17 12 13 os.chdir(buildroot) 18 14 19 os.chdir(buildroot)20 #execfile('test_all.py')21 15 -
anuga_work/development/anuga_1d/dam_h_elevation.py
r5737 r5738 42 42 43 43 domain.default_order = 2 44 domain.default_time_order = 244 domain.default_time_order = 1 45 45 domain.cfl = 1.0 46 46 domain.limiter = "vanleer" … … 61 61 for i in range(len(x)): 62 62 if 1100.0 <= x[i] <=1200.0: 63 y[i]=10.0 0163 y[i]=10.0 64 64 else: 65 65 y[i]=10.0 … … 140 140 from pylab import clf,plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ion 141 141 #print 'Test1' 142 ion() 142 143 hold(False) 143 144 #clf() … … 150 151 151 152 plot(X,ElevationQ,X,HeightQ) 152 plot1.set_ylim([ -1,12])153 plot1.set_ylim([9.99,10.01]) 153 154 xlabel('Position') 154 155 ylabel('Stage') … … 161 162 plot2 = subplot(212) 162 163 plot(X,MomentumQ) 163 plot2.set_ylim([-1,1])164 #plot2.set_ylim([-1,1]) 164 165 165 166 #legend( ('Numerical Solution', 'for momentum'), … … 173 174 file += ".eps" 174 175 #savefig(file) 175 show()176 #show() 176 177 177 178 print 'That took %.2f seconds'%(time.time()-t0) 179 180 raw_input("Press any key") -
anuga_work/development/anuga_1d/domain.py
r5737 r5738 697 697 self.datadir = name 698 698 699 #Main components of evolve 700 def evolve(self, yieldstep = None, finaltime = None, 699 700 #-------------------------- 701 # Main components of evolve 702 #-------------------------- 703 704 def evolve(self, 705 yieldstep = None, 706 finaltime = None, 707 duration = None, 708 skip_initial_step = False): 709 """Evolve model through time starting from self.starttime. 710 711 712 yieldstep: Interval between yields where results are stored, 713 statistics written and domain inspected or 714 possibly modified. If omitted the internal predefined 715 max timestep is used. 716 Internally, smaller timesteps may be taken. 717 718 duration: Duration of simulation 719 720 finaltime: Time where simulation should end. This is currently 721 relative time. So it's the same as duration. 722 723 If both duration and finaltime are given an exception is thrown. 724 725 726 skip_initial_step: Boolean flag that decides whether the first 727 yield step is skipped or not. This is useful for example to avoid 728 duplicate steps when multiple evolve processes are dove tailed. 729 730 731 Evolve is implemented as a generator and is to be called as such, e.g. 732 733 for t in domain.evolve(yieldstep, finaltime): 734 <Do something with domain and t> 735 736 737 All times are given in seconds 738 739 """ 740 741 from config import min_timestep, max_timestep, epsilon 742 743 # FIXME: Maybe lump into a larger check prior to evolving 744 msg = 'Boundary tags must be bound to boundary objects before ' 745 msg += 'evolving system, ' 746 msg += 'e.g. using the method set_boundary.\n' 747 msg += 'This system has the boundary tags %s '\ 748 %self.get_boundary_tags() 749 assert hasattr(self, 'boundary_objects'), msg 750 751 752 if yieldstep is None: 753 yieldstep = max_timestep 754 else: 755 yieldstep = float(yieldstep) 756 757 self._order_ = self.default_order 758 759 760 if finaltime is not None and duration is not None: 761 # print 'F', finaltime, duration 762 msg = 'Only one of finaltime and duration may be specified' 763 raise msg 764 else: 765 if finaltime is not None: 766 self.finaltime = float(finaltime) 767 if duration is not None: 768 self.finaltime = self.starttime + float(duration) 769 770 771 772 N = len(self) # Number of triangles 773 self.yieldtime = 0.0 # Track time between 'yields' 774 775 # Initialise interval of timestep sizes (for reporting only) 776 self.min_timestep = max_timestep 777 self.max_timestep = min_timestep 778 self.number_of_steps = 0 779 self.number_of_first_order_steps = 0 780 781 782 # Update ghosts 783 self.update_ghosts() 784 785 # Initial update of vertex and edge values 786 self.distribute_to_vertices_and_edges() 787 788 # Update extrema if necessary (for reporting) 789 self.update_extrema() 790 791 # Initial update boundary values 792 self.update_boundary() 793 794 # Or maybe restore from latest checkpoint 795 if self.checkpoint is True: 796 self.goto_latest_checkpoint() 797 798 if skip_initial_step is False: 799 yield(self.time) # Yield initial values 800 801 while True: 802 803 # Evolve One Step, using appropriate timestepping method 804 if self.get_timestepping_method() == 'euler': 805 self.evolve_one_euler_step(yieldstep,finaltime) 806 807 elif self.get_timestepping_method() == 'rk2': 808 self.evolve_one_rk2_step(yieldstep,finaltime) 809 810 elif self.get_timestepping_method() == 'rk3': 811 self.evolve_one_rk3_step(yieldstep,finaltime) 812 813 # Update extrema if necessary (for reporting) 814 self.update_extrema() 815 816 817 self.yieldtime += self.timestep 818 self.number_of_steps += 1 819 if self._order_ == 1: 820 self.number_of_first_order_steps += 1 821 822 # Yield results 823 if finaltime is not None and self.time >= finaltime-epsilon: 824 825 if self.time > finaltime: 826 # FIXME (Ole, 30 April 2006): Do we need this check? 827 # Probably not (Ole, 18 September 2008). Now changed to 828 # Exception 829 msg = 'WARNING (domain.py): time overshot finaltime. ' 830 msg += 'Contact Ole.Nielsen@ga.gov.au' 831 raise Exception, msg 832 833 834 # Yield final time and stop 835 self.time = finaltime 836 yield(self.time) 837 break 838 839 840 if self.yieldtime >= yieldstep: 841 # Yield (intermediate) time and allow inspection of domain 842 843 if self.checkpoint is True: 844 self.store_checkpoint() 845 self.delete_old_checkpoints() 846 847 # Pass control on to outer loop for more specific actions 848 yield(self.time) 849 850 # Reinitialise 851 self.yieldtime = 0.0 852 self.min_timestep = max_timestep 853 self.max_timestep = min_timestep 854 self.number_of_steps = 0 855 self.number_of_first_order_steps = 0 856 self.max_speed = zeros(N, Float) 857 858 859 def evolve_one_euler_step(self, yieldstep, finaltime): 860 """ 861 One Euler Time Step 862 Q^{n+1} = E(h) Q^n 863 """ 864 865 # Compute fluxes across each element edge 866 self.compute_fluxes() 867 868 # Update timestep to fit yieldstep and finaltime 869 self.update_timestep(yieldstep, finaltime) 870 871 # Update conserved quantities 872 self.update_conserved_quantities() 873 874 # Update ghosts 875 self.update_ghosts() 876 877 # Update vertex and edge values 878 self.distribute_to_vertices_and_edges() 879 880 # Update boundary values 881 self.update_boundary() 882 883 # Update time 884 self.time += self.timestep 885 886 887 888 889 def evolve_one_rk2_step(self, yieldstep, finaltime): 890 """ 891 One 2nd order RK timestep 892 Q^{n+1} = 0.5 Q^n + 0.5 E(h)^2 Q^n 893 """ 894 895 # Save initial initial conserved quantities values 896 self.backup_conserved_quantities() 897 898 #-------------------------------------- 899 # First euler step 900 #-------------------------------------- 901 902 # Compute fluxes across each element edge 903 self.compute_fluxes() 904 905 # Update timestep to fit yieldstep and finaltime 906 self.update_timestep(yieldstep, finaltime) 907 908 # Update conserved quantities 909 self.update_conserved_quantities() 910 911 # Update ghosts 912 self.update_ghosts() 913 914 # Update vertex and edge values 915 self.distribute_to_vertices_and_edges() 916 917 # Update boundary values 918 self.update_boundary() 919 920 # Update time 921 self.time += self.timestep 922 923 #------------------------------------ 924 # Second Euler step 925 #------------------------------------ 926 927 # Compute fluxes across each element edge 928 self.compute_fluxes() 929 930 # Update conserved quantities 931 self.update_conserved_quantities() 932 933 #------------------------------------ 934 # Combine initial and final values 935 # of conserved quantities and cleanup 936 #------------------------------------ 937 938 # Combine steps 939 self.saxpy_conserved_quantities(0.5, 0.5) 940 941 #----------------------------------- 942 # clean up vertex values 943 #----------------------------------- 944 945 # Update ghosts 946 self.update_ghosts() 947 948 # Update vertex and edge values 949 self.distribute_to_vertices_and_edges() 950 951 # Update boundary values 952 self.update_boundary() 953 954 955 956 def evolve_one_rk3_step(self, yieldstep, finaltime): 957 """ 958 One 3rd order RK timestep 959 Q^(1) = 3/4 Q^n + 1/4 E(h)^2 Q^n (at time t^n + h/2) 960 Q^{n+1} = 1/3 Q^n + 2/3 E(h) Q^(1) (at time t^{n+1}) 961 """ 962 963 # Save initial initial conserved quantities values 964 self.backup_conserved_quantities() 965 966 initial_time = self.time 967 968 #-------------------------------------- 969 # First euler step 970 #-------------------------------------- 971 972 # Compute fluxes across each element edge 973 self.compute_fluxes() 974 975 # Update timestep to fit yieldstep and finaltime 976 self.update_timestep(yieldstep, finaltime) 977 978 # Update conserved quantities 979 self.update_conserved_quantities() 980 981 # Update ghosts 982 self.update_ghosts() 983 984 # Update vertex and edge values 985 self.distribute_to_vertices_and_edges() 986 987 # Update boundary values 988 self.update_boundary() 989 990 # Update time 991 self.time += self.timestep 992 993 #------------------------------------ 994 # Second Euler step 995 #------------------------------------ 996 997 # Compute fluxes across each element edge 998 self.compute_fluxes() 999 1000 # Update conserved quantities 1001 self.update_conserved_quantities() 1002 1003 #------------------------------------ 1004 #Combine steps to obtain intermediate 1005 #solution at time t^n + 0.5 h 1006 #------------------------------------ 1007 1008 # Combine steps 1009 self.saxpy_conserved_quantities(0.25, 0.75) 1010 1011 # Update ghosts 1012 self.update_ghosts() 1013 1014 # Update vertex and edge values 1015 self.distribute_to_vertices_and_edges() 1016 1017 # Update boundary values 1018 self.update_boundary() 1019 1020 # Set substep time 1021 self.time = initial_time + self.timestep*0.5 1022 1023 #------------------------------------ 1024 # Third Euler step 1025 #------------------------------------ 1026 1027 # Compute fluxes across each element edge 1028 self.compute_fluxes() 1029 1030 # Update conserved quantities 1031 self.update_conserved_quantities() 1032 1033 #------------------------------------ 1034 # Combine final and initial values 1035 # and cleanup 1036 #------------------------------------ 1037 1038 # Combine steps 1039 self.saxpy_conserved_quantities(2.0/3.0, 1.0/3.0) 1040 1041 # Update ghosts 1042 self.update_ghosts() 1043 1044 # Update vertex and edge values 1045 self.distribute_to_vertices_and_edges() 1046 1047 # Update boundary values 1048 self.update_boundary() 1049 1050 # Set new time 1051 self.time = initial_time + self.timestep 1052 1053 1054 def backup_conserved_quantities(self): 1055 N = len(self) # Number_of_triangles 1056 1057 # Backup conserved_quantities centroid values 1058 for name in self.conserved_quantities: 1059 Q = self.quantities[name] 1060 Q.backup_centroid_values() 1061 1062 def saxpy_conserved_quantities(self,a,b): 1063 N = len(self) #number_of_triangles 1064 1065 # Backup conserved_quantities centroid values 1066 for name in self.conserved_quantities: 1067 Q = self.quantities[name] 1068 Q.saxpy_centroid_values(a,b) 1069 1070 1071 #============================== 1072 # John Jakeman's old evolve code 1073 #============================= 1074 1075 def evolve_john(self, yieldstep = None, finaltime = None, 701 1076 skip_initial_step = False): 702 1077 """Evolve model from time=0.0 to finaltime yielding results … … 728 1103 729 1104 self.order = self.default_order 1105 730 1106 self.time_order = self.default_time_order 731 1107 … … 741 1117 742 1118 #update ghosts 743 #self.update_ghosts()1119 self.update_ghosts() 744 1120 745 1121 #Initial update of vertex and edge values … … 969 1345 970 1346 1347 def update_ghosts(self): 1348 pass 1349 971 1350 def update_boundary(self): 972 1351 """Go through list of boundary objects and update boundary values … … 1035 1414 1036 1415 self.timestep = timestep 1037 1416 1417 def update_extrema(self): 1418 pass 1038 1419 1039 1420 def compute_forcing_terms(self): -
anuga_work/development/anuga_1d/quantity.py
r5536 r5738 51 51 #Allocate space for other quantities 52 52 self.centroid_values = zeros(N, Float) 53 self.centroid_backup_values = zeros(N, Float) 53 54 #self.edge_values = zeros((N, 2), Float) 54 55 #edge values are values of the ends of each interval … … 854 855 else: 855 856 qv[k,i] = qc[k] 856 857 858 def backup_centroid_values(self): 859 # Call correct module function 860 # (either from this module or C-extension) 861 #backup_centroid_values(self) 862 N = self.domain.number_of_elements 863 for i in range(self.N): 864 quantity.centroid_backup_values[i] = quantity.centroid_values[i] 865 866 def saxpy_centroid_values(self,a,b): 867 # Call correct module function 868 # (either from this module or C-extension) 869 #saxpy_centroid_values(self,a,b) 870 N = self.domain.number_of_elements 871 for i in range(self.N): 872 quantity.centroid_backup_values[i] = a*quantity.centroid_values[i] + b*quantity.centroid_backup_values[i] 857 873 858 874 class Conserved_quantity(Quantity): … … 869 885 print "Use Quantity instead of Conserved_quantity" 870 886 871 887 872 888 873 889 def newLinePlot(title='Simple Plot'): -
anuga_work/development/anuga_1d/test_shallow_water.py
r5731 r5738 41 41 #print domain.quantities['stage'].explicit_update 42 42 #print domain.quantities['xmomentum'].explicit_update 43 print stage_ud43 #print stage_ud 44 44 45 45 assert allclose( domain.quantities['stage'].explicit_update, stage_ud ) … … 134 134 domain.default_order = 1 135 135 domain.default_time_order = 1 136 yieldstep=1 0.0137 finaltime=1 0.0136 yieldstep=1.0 137 finaltime=1.0 138 138 139 139 for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime): … … 168 168 domain.default_order = 2 169 169 domain.default_time_order = 1 170 yieldstep=1 0.0171 finaltime=1 0.0170 yieldstep=1.0 171 finaltime=1.0 172 172 173 173 for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime): … … 177 177 assert allclose( zeros(10), domain.quantities['xmomentum'].centroid_values ) 178 178 179 def test_evolve_ euler_second_order_space_time(self):179 def test_evolve_second_order_space_time(self): 180 180 """ 181 181 Compare still lake solution for various versions of shallow_water_domain … … 192 192 domain.default_order = 2 193 193 domain.default_time_order = 2 194 yieldstep=1 0.0195 finaltime=1 0.0194 yieldstep=1.0 195 finaltime=1.0 196 196 197 197 for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
Note: See TracChangeset
for help on using the changeset viewer.