Last change on this file since 5533 was 4959, checked in by steve, 16 years ago

Copied John Jakeman's and Sudi Mungkasi's 1d shallow water code

File size: 6.9 KB
Line
1import os
2from math import sqrt, pi, sin, cos
4from Numeric import allclose, array, zeros, ones, Float, take, sqrt
5from config import g, epsilon
6
7def newLinePlot(title='Simple Plot'):
8    import Gnuplot
9    g = Gnuplot.Gnuplot(persist=0)
10    g.title(title)
11    g('set data style linespoints')
12    g.xlabel('x')
13    g.ylabel('y')
14    return g
15
16def linePlot(g,x1,y1,x2,y2,x3,y3):
17    import Gnuplot
18    plot1 = Gnuplot.PlotItems.Data(x1.flat,y1.flat,with="lines 2")
19    plot2 = Gnuplot.PlotItems.Data(x2.flat,y2.flat,with="lines 3")
20    plot3 = Gnuplot.PlotItems.Data(x3.flat,y3.flat,with="linespoints 1")
21    g.plot(plot1,plot2,plot3)
22
23
24def analytic_cannal(C,t):
25    N = len(C)
26    u = zeros(N,Float)    ## water velocity
27    h = zeros(N,Float)    ## water depth
28    x = C
29    g = 9.81
30
31
32    ## Define Basin Bathymetry
33    z_b = zeros(N,Float) ## elevation of basin
34    z = zeros(N,Float)   ## elevation of water surface
35    z_infty = 10.0       ## max equilibrium water depth at lowest point.
36    L_x = 2500.0         ## width of channel
37
38    A0 = 0.5*L_x                  ## determines amplitudes of oscillations
39    omega = sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation
40
41    #x1 = A0*cos(omega*t)-L_x # left shoreline
42    #x2 = A0*cos(omega*t)+L_x # right shoreline
43
44    for i in range(N):
45        z_b[i] = z_infty*(x[i]**2/L_x**2)
46        u[i] = -A0*omega*sin(omega*t)
47        z[i] = z_infty+2*A0*z_infty/L_x*cos(omega*t)*(x[i]/L_x-0.5*A0/(L_x)*cos(omega*t))
48
49    h = z-z_b
50    return u,h,z,z_b
51
52
53#plot2 = newLinePlot("Momentum")
54
55L_x = 2500.0     # Length of channel (m)
56N = 400         # Number of compuational cells
57cell_len = 4*L_x/N   # Origin = 0.0
58
59points = zeros(N+1,Float)
60for i in range(N+1):
61        points[i] = -2*L_x +i*cell_len
62
63domain = Domain(points)
64
65domain.order = 2 #make this unnecessary
66domain.default_order = 2
67domain.default_time_order = 1
68domain.cfl = 1.0
69domain.beta = 1.0
70#domain.limiter = "minmod_kurganov"
71domain.limiter = "vanleer"
73#domain.limiter = "superbee"
74#domain.limiter = "minmod"
75
76def stage(x):
77
78    z_infty = 10.0       ## max equilibrium water depth at lowest point.
79    L_x = 2500.0         ## width of channel
80    A0 = 0.5*L_x                  ## determines amplitudes of oscillations
81    omega = sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation
82    t=0.0
83    y = zeros(len(x),Float)
84    for i in range(len(x)):
85        y[i] = z_infty+2*A0*z_infty/L_x*cos(omega*t)*(x[i]/L_x-0.5*A0/(L_x)*cos(omega*t))
86    return y
87
88def elevation(x):
89    N = len(x)
90    z_infty = 10.0
91    z = zeros(N,Float)
92    L_x = 2500.0
93    A0 = 0.5*L_x
94    omega = sqrt(2*g*z_infty)/L_x
95    for i in range(N):
96        z[i] = z_infty*(x[i]**2/L_x**2)
97    return z
98
99def height(x):
100    z_infty = 10.0       ## max equilibrium water depth at lowest point.
101    L_x = 2500.0         ## width of channel
102    A0 = 0.5*L_x                  ## determines amplitudes of oscillations
103    omega = sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation
104    t=0.0
105    y = zeros(len(x),Float)
106    for i in range(len(x)):
107        y[i] = max(z_infty+2*A0*z_infty/L_x*cos(omega*t)*(x[i]/L_x-0.5*A0/(L_x)*cos(omega*t))-z_infty*(x[i]**2/L_x**2),0.0)
108    return y
109
110
111domain.set_quantity('stage', stage)
112domain.set_quantity('elevation',elevation)
113#domain.set_quantity('height',height)
114
115domain.set_boundary({'exterior': Reflective_boundary(domain)})
116
117import time
118t0 = time.time()
119#yieldstep = 50.0
120#yieldstep = 10.0
121#yieldstep = 1122.0
122
123#finaltime = 1122.0*0.999488
124finaltime=0.75*1122.0
125yieldstep = finaltime
126#yieldstep = finaltime
127#finaltime = 1122.0*4.0
128#finaltime = 2000.0
129
130plot1 = newLinePlot("Stage")
131plot2 = newLinePlot("Xmomentum")
132C = domain.centroids
133X = domain.vertices
134
135ElevationQ = domain.quantities['elevation'].vertex_values
136StageQ = domain.quantities['stage'].vertex_values
137XmomQ = domain.quantities['xmomentum'].vertex_values
138#HeightC = domain.quantities['stage'].centroid_values-domain.quantities['elevation'].centroid_values
139HeightC = domain.quantities['height'].centroid_values
140XmomC = domain.quantities['xmomentum'].centroid_values
141from util import calculate_wetted_area
142import time
143
144t0 = time.time()
145for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
146    pass
147    #print StageQ
148    #print "integral %f"%domain.quantities['height'].get_integral()
149    #if t == 0.0:
150    #     initial_integral = domain.quantities['height'].get_integral()
151    #assert(allclose(initial_integral,domain.quantities['height'].get_integral()))
152    #domain.write_time()
153    u,h,z,z_b = analytic_cannal(X.flat,t)
154    linePlot(plot1,X,z,X,z_b,X,StageQ)
155    linePlot(plot2,X,u*h,X,u*h,X,XmomQ)
156    #HeightQ = domain.quantities['stage'].centroid_values-domain.quantities['elevation'].centroid_values
157    #u,hc,z,z_b = analytic_cannal(C,t)
158    #for k in range(N):
159    #    if hc[k] < 0.0:
160    #        hc[k] = 0.0
161    #error = 1.0/(N)*sum(abs(hc-HeightQ))
162    #print 'Error measured at centroids', error
163
164#assert(allclose(initial_integral,domain.quantities['height'].get_integral()))
165
166u,h,z,z_b = analytic_cannal(C.flat,t)
167z_infty = 10.0
168L_x = 2500.0
169A0 = 0.5*L_x
170omega = sqrt(2*g*z_infty)/L_x
171x1 = A0*cos(omega*t)-L_x # left shoreline
172x2 = A0*cos(omega*t)+L_x # right shoreline
173for j in range(N):
174    if (C[j]<x1) | (C[j] > x2):
175            u[j] = 0.0
176            h[j] = 0.0
177#print h
178#print HeightC
179h_error = 1.0/(N)*sum(abs(h-HeightC))
180uh_error = 1.0/(N)*sum(abs(u*h-XmomC))
181print "h_error %.10f" %(h_error)
182print "uh_error %.10f"% (uh_error)
183
184from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,rc
185rc('text', usetex=True)
186X = domain.vertices
187u,h,z,z_b = analytic_cannal(X.flat,domain.time)
188StageQ = domain.quantities['stage'].vertex_values
189XmomQ = domain.quantities['xmomentum'].vertex_values
190hold(False)
191plot1 = subplot(211)
192z_infty = 10.0
193L_x = 2500.0
194A0 = 0.5*L_x
195omega = sqrt(2*g*z_infty)/L_x
196x1 = A0*cos(omega*t)-L_x # left shoreline
197x2 = A0*cos(omega*t)+L_x # right shoreline
198for n in range(N):
199    m = 2*n
200    for j in range(2):
201        if (X.flat[m+j]<x1) | (X.flat[m+j] > x2):
202            u[m+j] = 0.0
203plot(X,z,X,StageQ,X,z_b)
204#plot1.set_ylim([4,11])
205#title('Free Surface Elevation of a Dry Dam-Break')
206#xlabel('x (m)')
207ylabel('Stage (m)')
208legend(('Analytical Solution', 'Numerical Solution'),
210plot2 = subplot(212)
211plot(X,u*h,X,XmomQ)
212#plot2.set_ylim([-1,25])
213#title('Xmomentum Profile of a Dry Dam-Break')
214xlabel('x (m)')
215ylabel(r'x-momentum (\$m^2/s\$)')
217filename += domain.limiter
218filename += str(finaltime)
219filename += ".eps"
220print filename
221savefig(filename)
222show()
223
224"""
226h_error 0.0620113496
227uh_error 0.7105544340