1 | import os |
---|
2 | import Gnuplot |
---|
3 | from math import sqrt, pi |
---|
4 | from shallow_water_1d import * |
---|
5 | from Numeric import allclose, array, zeros, ones, Float, take, sqrt |
---|
6 | from config import g, epsilon |
---|
7 | |
---|
8 | |
---|
9 | #h_0 = 0.48 |
---|
10 | #u_0 = 0.625 |
---|
11 | #h_1 = 0.48 |
---|
12 | #u_1 = 0.625 |
---|
13 | |
---|
14 | #h_0 = 0.48 |
---|
15 | #u_0 = 0.625 |
---|
16 | #h_1 = 0.106 |
---|
17 | #u_1 = 2.829228 |
---|
18 | |
---|
19 | h_0 = 0.5 |
---|
20 | u_0 = 0.6 |
---|
21 | h_1 = 0.5 |
---|
22 | u_1 = 0.6 |
---|
23 | g = 9.81 |
---|
24 | |
---|
25 | #h_c = (q**2/g)**(1/3) |
---|
26 | #u_c = sqrt(g*h_c) |
---|
27 | |
---|
28 | q = 0.3 |
---|
29 | g = 9.81 |
---|
30 | |
---|
31 | b = 4.0 |
---|
32 | xmax = 25.0 |
---|
33 | z_bmax = 0.2 |
---|
34 | Fr_0 = u_0/sqrt(g*h_0) |
---|
35 | |
---|
36 | def stage(x): |
---|
37 | y,u,h = analytical_sol(x) |
---|
38 | return y |
---|
39 | |
---|
40 | |
---|
41 | def elevation(x): |
---|
42 | b = 4.0 |
---|
43 | xmax = 25.0 |
---|
44 | z_bmax = 0.2 |
---|
45 | z_b = zeros(len(x),Float) |
---|
46 | for i in range(len(x)): |
---|
47 | if (x[i] >= xmax-b) & (x[i] <= xmax+b): |
---|
48 | z_b[i] = z_bmax*(1.0-(x[i]-xmax)**2/b**2) |
---|
49 | else: |
---|
50 | z_b[i] = 0.0 |
---|
51 | return z_b |
---|
52 | #return 0.0 |
---|
53 | |
---|
54 | def xmomentum(x): |
---|
55 | return u_0*h_0 |
---|
56 | |
---|
57 | |
---|
58 | def bisection(func,xR,xL,H): |
---|
59 | while ((xR - xL) > epsilon): |
---|
60 | |
---|
61 | #Calculate midpoint of domain |
---|
62 | xM = xL + (xR - xL) / 2.0 |
---|
63 | |
---|
64 | #Find f(xM) |
---|
65 | if ((f(xL,H) * f(xM,H)) > 0): |
---|
66 | #Throw away left half |
---|
67 | xL = xM |
---|
68 | else: |
---|
69 | #Throw away right half |
---|
70 | xR = xM |
---|
71 | return xR |
---|
72 | |
---|
73 | def f(D,H): |
---|
74 | return D**3+D**2*(H-1.0-Fr_0**2/2.0)+Fr_0**2/2.0 |
---|
75 | |
---|
76 | def fprime(D,H): |
---|
77 | Fr_0 = u_0/sqrt(g*h_0) |
---|
78 | return 3*D**2+2*D*(H-1-Fr_0**2/2.0) |
---|
79 | |
---|
80 | def analytical_sol(x): |
---|
81 | y = zeros(len(x),Float) |
---|
82 | u = zeros(len(x),Float) |
---|
83 | height = zeros(len(x),Float) |
---|
84 | for i in range(len(x)): |
---|
85 | if (x[i] >= xmax-b) & (x[i] <= xmax+b): |
---|
86 | zb = z_bmax*(1.0-(x[i]-xmax)**2/b**2) |
---|
87 | else: |
---|
88 | zb = 0.0 |
---|
89 | |
---|
90 | H = zb/h_0 |
---|
91 | y1 = bisection(f,1.0,0.2/0.5,H) |
---|
92 | height[i] = y1*h_0 |
---|
93 | y[i] = height[i]+zb |
---|
94 | u[i] = q/(height[i]) |
---|
95 | |
---|
96 | |
---|
97 | return y,u,height |
---|
98 | |
---|
99 | Fr_0 = u_0/sqrt(g*h_0) |
---|
100 | zb=0.2 |
---|
101 | H = zb/h_0 |
---|
102 | D = bisection(f,1,0.2/0.5,H) |
---|
103 | |
---|
104 | L = 50.0 |
---|
105 | N = 50 #800 |
---|
106 | cell_len = L/N |
---|
107 | points = zeros(N+1,Float) |
---|
108 | for j in range(N+1): |
---|
109 | points[j] = j*cell_len |
---|
110 | boundary = { (0,0): 'left',(N-1,1): 'right'} |
---|
111 | domain = Domain(points,boundary) |
---|
112 | D1 = Dirichlet_boundary([h_0,u_0*h_0]) |
---|
113 | D2 = Dirichlet_boundary([h_1,u_1*h_1]) |
---|
114 | domain.set_boundary({'left':D1,'right':D2}) |
---|
115 | domain.set_quantity('elevation',elevation) |
---|
116 | domain.set_quantity('stage',stage) |
---|
117 | domain.set_quantity('xmomentum',xmomentum) |
---|
118 | X = domain.vertices |
---|
119 | C = domain.centroids |
---|
120 | |
---|
121 | Stage = domain.quantities['stage'] |
---|
122 | Xmom = domain.quantities['xmomentum'] |
---|
123 | |
---|
124 | StageQ = Stage.vertex_values |
---|
125 | StageC = Stage.centroid_values |
---|
126 | XmomQ = Xmom.vertex_values |
---|
127 | XmomC = Xmom.centroid_values |
---|
128 | ElevationQ = domain.quantities['elevation'].vertex_values |
---|
129 | |
---|
130 | stage_bdry = Stage.boundary_values |
---|
131 | xmom_bdry = Xmom.boundary_values |
---|
132 | |
---|
133 | w,u,h = analytical_sol(X.flat) |
---|
134 | wc, uc, hc = analytical_sol(C) |
---|
135 | |
---|
136 | from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ylim,xlim, rc |
---|
137 | hold(False) |
---|
138 | rc('text', usetex=True) |
---|
139 | h_error = zeros(1,Float) |
---|
140 | uh_error = zeros(1,Float) |
---|
141 | k = 0 |
---|
142 | yieldstep = 25.0 |
---|
143 | finaltime = 25.0 |
---|
144 | domain.limiter = "vanleer" |
---|
145 | #domain.limiter = "steve_minmod" |
---|
146 | domain.default_order = 2 |
---|
147 | domain.default_time_order = 2#check order is not working |
---|
148 | for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): |
---|
149 | domain.write_time() |
---|
150 | print 'Test 1' |
---|
151 | plot1 = subplot(211) |
---|
152 | print 'Test 2' |
---|
153 | plot(X,w,X,StageQ,X,ElevationQ) |
---|
154 | print 'Test 3' |
---|
155 | #xlabel('Position') |
---|
156 | #ylabel('Stage (m)') |
---|
157 | #legend(('Analytical Solution', 'Numerical Solution', 'Channel Bed'), |
---|
158 | # 'upper right', shadow=True) |
---|
159 | plot1.set_ylim([-0.1,1.0]) |
---|
160 | print 'Test 4' |
---|
161 | plot2 = subplot(212) |
---|
162 | print 'Test 5' |
---|
163 | plot(X,u*h,X,XmomQ)#/(StageQ-ElevationQ)) |
---|
164 | print 'Test 6' |
---|
165 | #xlabel('x (m)') |
---|
166 | #ylabel(r'X-momentum ($m^2/s$)') |
---|
167 | h_error[k] = 1.0/(N)*sum(abs(wc-StageC)) |
---|
168 | uh_error[k] = 1.0/(N)*sum(abs(uc*hc-XmomC)) |
---|
169 | print "h_error %.10f" %(h_error[k]) |
---|
170 | print "uh_error %.10f"% (uh_error[k]) |
---|
171 | print "h_max %.10f"%max(wc-StageC) |
---|
172 | print "uh max %.10f"%max(uc*hc-XmomC) |
---|
173 | filename = "steady_flow_" |
---|
174 | filename += domain.limiter |
---|
175 | filename += str(400) |
---|
176 | filename += ".ps" |
---|
177 | #savefig(filename) |
---|
178 | show() |
---|
179 | |
---|