 Timestamp:
 Jun 16, 2010, 5:30:17 PM (13 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_work/development/2010projects/anuga_1d/base/quantity.py
r7850 r7852 838 838 second order scheme. 839 839 """ 840 840 841 if self.domain.limiter == "pyvolution": 841 #Z = self.gradients842 #print "gradients 1",Z843 self.compute_gradients()844 #print "gradients 2",Z845 846 #Z = self.gradients847 #print "gradients 1",Z848 #self.compute_minmod_gradients()849 #print "gradients 2", Z850 851 G = self.gradients852 V = self.domain.vertices853 qc = self.centroid_values854 qv = self.vertex_values855 856 #Check each triangle857 for k in range(self.domain.number_of_elements):858 #Centroid coordinates859 x = self.domain.centroids[k]860 861 #vertex coordinates862 x0, x1 = V[k,:]863 864 #Extrapolate865 qv[k,0] = qc[k] + G[k]*(x0x)866 qv[k,1] = qc[k] + G[k]*(x1x)867 842 self.limit_pyvolution() 868 843 869 844 elif self.domain.limiter == "minmod_steve": 870 845 self.limit_minmod() 846 871 847 else: 872 848 self.limit_range() 873 849 874 875 876 def limit_minmod(self): 877 #Z = self.gradients 878 #print "gradients 1",Z 879 self.compute_minmod_gradients() 880 #print "gradients 2", Z 881 882 G = self.gradients 883 V = self.domain.vertices 884 qc = self.centroid_values 885 qv = self.vertex_values 886 887 #Check each triangle 888 for k in range(self.domain.number_of_elements): 889 #Centroid coordinates 890 x = self.domain.centroids[k] 891 892 #vertex coordinates 893 x0, x1 = V[k,:] 894 895 #Extrapolate 896 qv[k,0] = qc[k] + G[k]*(x0x) 897 qv[k,1] = qc[k] + G[k]*(x1x) 898 899 900 def limit_pyvolution(self): 901 """ 902 Limit slopes for each volume to eliminate artificial variance 903 introduced by e.g. second order extrapolator 904 905 This is an unsophisticated limiter as it does not take into 906 account dependencies among quantities. 907 908 precondition: 909 vertex values are estimated from gradient 910 postcondition: 911 vertex values are updated 912 """ 913 914 915 N = self.domain.number_of_elements 916 beta = self.domain.beta 917 #beta = 0.8 918 919 qc = self.centroid_values 920 qv = self.vertex_values 921 922 #Find min and max of this and neighbour's centroid values 850 851 852 853 854 def find_qmax_qmin(self): 855 """ Find min and max of this and neighbour's centroid values""" 856 857 923 858 qmax = self.qmax 924 859 qmin = self.qmin 860 861 qc = self.centroid_values 925 862 926 863 for k in range(N): … … 935 872 936 873 874 875 def limit_minmod(self): 876 #Z = self.gradients 877 #print "gradients 1",Z 878 self.compute_minmod_gradients() 879 #print "gradients 2", Z 880 881 G = self.gradients 882 V = self.domain.vertices 883 qc = self.centroid_values 884 qv = self.vertex_values 885 886 x = self.domain.centroids 887 888 x0 = V[:,0] 889 x1 = V[:,1] 890 891 #Extrapolate 892 qv[:,0] = qc + G*(x0x) 893 qv[:,1] = qc + G*(x1x) 894 895 # #Check each triangle 896 # for k in range(self.domain.number_of_elements): 897 # #Centroid coordinates 898 # x = self.domain.centroids[k] 899 # 900 # #vertex coordinates 901 # x0, x1 = V[k,:] 902 # 903 # #Extrapolate 904 # qv[k,0] = qc[k] + G[k]*(x0x) 905 # qv[k,1] = qc[k] + G[k]*(x1x) 906 907 908 def limit_pyvolution(self): 909 """ 910 Limit slopes for each volume to eliminate artificial variance 911 introduced by e.g. second order extrapolator 912 913 This is an unsophisticated limiter as it does not take into 914 account dependencies among quantities. 915 916 precondition: 917 vertex values are estimated from gradient 918 postcondition: 919 vertex values are updated 920 """ 921 922 923 N = self.domain.number_of_elements 924 beta = self.domain.beta 925 #beta = 0.8 926 927 self.compute_gradients() 928 929 930 G = self.gradients 931 V = self.domain.vertices 932 C = self.domain.centroids 933 qc = self.centroid_values 934 qv = self.vertex_values 935 936 V0 = V[:,0] 937 V1 = V[:,1] 938 939 # Extrapolate to vertices 940 qv[:,0] = qc + G*(V0C) 941 qv[:,1] = qc + G*(V1C) 942 943 944 # Find max and min values 945 self.find_qmax_qmin() 946 947 qmax = self.qmax 948 qmin = self.qmin 949 937 950 #Diffences between centroids and maxima/minima 938 951 dqmax = qmax  qc … … 941 954 #Deltas between vertex and centroid values 942 955 dq = numpy.zeros(qv.shape, numpy.float) 943 for i in range(2): 944 dq[:,i] = qv[:,i]  qc 945 946 #Phi limiter 947 for k in range(N): 948 949 #Find the gradient limiter (phi) across vertices 950 phi = 1.0 951 for i in range(2): 952 r = 1.0 953 if (dq[k,i] > 0): r = dqmax[k]/dq[k,i] 954 if (dq[k,i] < 0): r = dqmin[k]/dq[k,i] 955 956 phi = min( min(r*beta, 1), phi ) 957 958 #Then update using phi limiter 959 for i in range(2): 960 qv[k,i] = qc[k] + phi*dq[k,i] 956 957 dq[:,0] = qv[:,0]  qc 958 dq[:,1] = qv[:,1]  qc 959 960 phi = numpy.ones(qc.shape, numpy.float) 961 962 r0 = numpy.where(dq[:,0]>0.0,dqmax/dq[:,0],1.0) 963 r0 = numpy.where(dq[:,0]<0.0,dqmin/dq[:,0],r0) 964 965 r1 = numpy.where(dq[:,1]>0.0,dqmax/dq[:,1],1.0) 966 r1 = numpy.where(dq[:,1]<0.0,dqmin/dq[:,1],r1) 967 968 phi = numpy.min(r0*beta,numpy.min(r1*beta,1.0)) 969 970 qv[:,0] = qc + phi*dq[:,0] 971 qv[:,1] = qc + phi*dq[:,1] 972 973 # #Phi limiter 974 # for k in range(N): 975 # 976 # #Find the gradient limiter (phi) across vertices 977 # phi = 1.0 978 # for i in range(2): 979 # r = 1.0 980 # if (dq[k,i] > 0): r = dqmax[k]/dq[k,i] 981 # if (dq[k,i] < 0): r = dqmin[k]/dq[k,i] 982 # 983 # phi = min( min(r*beta, 1), phi ) 984 # 985 # #Then update using phi limiter 986 # for i in range(2): 987 # qv[k,i] = qc[k] + phi*dq[k,i] 961 988 962 989 def limit_range(self): 963 990 import sys 964 991 965 from utilimport minmod, minmod_kurganov, maxmod, vanleer, vanalbada992 from limiters_python import minmod, minmod_kurganov, maxmod, vanleer, vanalbada 966 993 967 994 limiter = self.domain.limiter … … 972 999 qc = self.centroid_values 973 1000 qv = self.vertex_values 974 C = self.domain.centroids 975 X = self.domain.vertices 1001 xc = self.domain.centroids 1002 x0 = self.domain.vertices[:,0] 1003 x1 = self.domain.vertices[:,1] 1004 976 1005 beta_p = numpy.zeros(N,numpy.float) 977 1006 beta_m = numpy.zeros(N,numpy.float) 978 1007 beta_x = numpy.zeros(N,numpy.float) 979 1008 980 for k in range(N): 981 982 n0 = self.domain.neighbours[k,0] 983 n1 = self.domain.neighbours[k,1] 1009 # for k in range(N): 1010 # 1011 # n0 = self.domain.neighbours[k,0] 1012 # n1 = self.domain.neighbours[k,1] 1013 # 1014 # if ( n0 >= 0) & (n1 >= 0): 1015 # #SLOPE DERIVATIVE LIMIT 1016 # beta_p[k] = (qc[k]qc[k1])/(C[k]C[k1]) 1017 # beta_m[k] = (qc[k+1]qc[k])/(C[k+1]C[k]) 1018 # beta_x[k] = (qc[k+1]qc[k1])/(C[k+1]C[k1]) 1019 1020 1021 1022 n0 = self.domain.neighbours[:,0] 1023 n1 = self.domain.neighbours[:,1] 1024 1025 1026 beta_p[1:] = (qc[1:]qc[:1])/(xc[1:]xc[:1]) 1027 beta_m[:1] = beta_p[1:] 1028 beta_x[1:1] = (qc[2:]qc[:2])/(xc[2:]xc[:2]) 1029 1030 dx = numpy.zeros(qv.shape, numpy.float) 1031 1032 dx[:,0] = x0  xc 1033 dx[:,1] = x1  xc 1034 1035 1036 phi = numpy.zeros(qc.shape, numpy.float) 1037 1038 phi[0] = (qc[1]  qc[0])/(xc[1]  xc[0]) 1039 phi[N1] = (qc[N1]  qc[N2])/(xc[N1]  xc[N2]) 1040 1041 if limiter == "minmod": 1042 phi[1:1] = minmod(beta_p[1:1],beta_m[1:1]) 1043 1044 elif limiter == "vanleer": 1045 phi[1:1] = vanleer(beta_p[1:1],beta_m[1:1]) 984 1046 985 if ( n0 >= 0) & (n1 >= 0): 986 #SLOPE DERIVATIVE LIMIT 987 beta_p[k] = (qc[k]qc[k1])/(C[k]C[k1]) 988 beta_m[k] = (qc[k+1]qc[k])/(C[k+1]C[k]) 989 beta_x[k] = (qc[k+1]qc[k1])/(C[k+1]C[k1]) 990 991 dq = numpy.zeros(qv.shape, numpy.float) 992 for i in range(2): 993 dq[:,i] =self.domain.vertices[:,i]self.domain.centroids 994 1047 elif limiter == "vanalbada": 1048 phi[1:1] = vanalbada(beta_p[1:1],beta_m[1:1]) 1049 1050 elif limiter == "minmod_kurganov": 1051 theta = 2.0 1052 phi[1:1] = minmod_kurganov(theta*beta_p[1:1],theta*beta_m[1:1], beta_x[1:1]) 1053 1054 elif limiter == "superbee": 1055 1056 slope1 = minmod(beta_m[1:1],2.0*beta_p[1:1]) 1057 slope2 = minmod(2.0*beta_m[1:1],beta_p[1:1]) 1058 phi[1:1] = maxmod(slope1,slope2) 1059 1060 1061 else: 1062 msg = 'Unknown limiter' 1063 raise Exception, msg 1064 1065 qv[:,0] = qc + phi*dx[:,0] 1066 qv[:,1] = qc + phi*dx[:,1] 1067 1068 995 1069 #Phi limiter 996 for k in range(N): 997 n0 = self.domain.neighbours[k,0] 998 n1 = self.domain.neighbours[k,1] 999 if n0 < 0: 1000 phi = (qc[k+1]  qc[k])/(C[k+1]  C[k]) 1001 elif n1 < 0: 1002 phi = (qc[k]  qc[k1])/(C[k]  C[k1]) 1003 #elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2): 1004 # phi = 0.0 1005 else: 1006 if limiter == "minmod": 1007 phi = minmod(beta_p[k],beta_m[k]) 1008 1009 elif limiter == "minmod_kurganov":#Change this 1010 # Also known as monotonized central difference limiter 1011 # if theta = 2.0 1012 theta = 2.0 1013 phi = minmod_kurganov(theta*beta_p[k],theta*beta_m[k],beta_x[k]) 1014 elif limiter == "superbee": 1015 slope1 = minmod(beta_m[k],2.0*beta_p[k]) 1016 slope2 = minmod(2.0*beta_m[k],beta_p[k]) 1017 phi = maxmod(slope1,slope2) 1018 1019 elif limiter == "vanleer": 1020 phi = vanleer(beta_p[k],beta_m[k]) 1021 1022 elif limiter == "vanalbada": 1023 phi = vanalbada(beta_m[k],beta_p[k]) 1024 1025 for i in range(2): 1026 qv[k,i] = qc[k] + phi*dq[k,i] 1070 # for k in range(N): 1071 # n0 = self.domain.neighbours[k,0] 1072 # n1 = self.domain.neighbours[k,1] 1073 # if n0 < 0: 1074 # phi = (qc[k+1]  qc[k])/(xc[k+1]  xc[k]) 1075 # elif n1 < 0: 1076 # phi = (qc[k]  qc[k1])/(xc[k]  xc[k1]) 1077 # #elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2): 1078 # # phi = 0.0 1079 # else: 1080 # if limiter == "minmod": 1081 # phi = minmod(beta_p[k],beta_m[k]) 1082 # 1083 # elif limiter == "minmod_kurganov":#Change this 1084 # # Also known as monotonized central difference limiter 1085 # # if theta = 2.0 1086 # theta = 2.0 1087 # phi = minmod_kurganov(theta*beta_p[k],theta*beta_m[k],beta_x[k]) 1088 # 1089 # elif limiter == "superbee": 1090 # slope1 = minmod(beta_m[k],2.0*beta_p[k]) 1091 # slope2 = minmod(2.0*beta_m[k],beta_p[k]) 1092 # phi = maxmod(slope1,slope2) 1093 # 1094 # elif limiter == "vanleer": 1095 # phi = vanleer(beta_p[k],beta_m[k]) 1096 # 1097 # elif limiter == "vanalbada": 1098 # phi = vanalbada(beta_m[k],beta_p[k]) 1099 # 1100 # for i in range(2): 1101 # qv[k,i] = qc[k] + phi*dx[k,i] 1027 1102 1028 1103 def limit_steve_slope(self): … … 1107 1182 1108 1183 1109 def newLinePlot(title='Simple Plot'): 1110 import pylab as g 1111 g.ion() 1112 g.hold(False) 1113 g.title(title) 1114 g.xlabel('x') 1115 g.ylabel('y') 1116 1117 1118 def linePlot(x,y): 1119 import pylab as g 1120 g.plot(x.flat,y.flat) 1121 1122 1123 def closePlots(): 1124 import pylab as g 1125 g.close('all') 1184 1126 1185 1127 1186 if __name__ == "__main__": 1128 1187 #from domain import Domain 1129 1188 from generic_domain import Generic_domain as Domain 1189 1190 1191 def newLinePlot(title='Simple Plot'): 1192 import pylab as g 1193 g.ion() 1194 g.hold(False) 1195 g.title(title) 1196 g.xlabel('x') 1197 g.ylabel('y') 1198 1199 1200 def linePlot(x,y): 1201 import pylab as g 1202 g.plot(x.flat,y.flat) 1203 1204 1205 def closePlots(): 1206 import pylab as g 1207 g.close('all') 1208 1209 1130 1210 1131 1211 points1 = [0.0, 1.0, 2.0, 3.0] … … 1213 1293 1214 1294 raw_input('press return to quit') 1215 closePlots() 1295 1296 closePlots()
Note: See TracChangeset
for help on using the changeset viewer.