source: anuga_work/development/anuga_1d/parabolic_cannal.py @ 5827

Last change on this file since 5827 was 5827, checked in by sudi, 15 years ago

Adding new versions of shallow_water_domain

File size: 5.1 KB
Line 
1import os
2from math import sqrt, pi, sin, cos
3from shallow_water_domain import *
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 = 100         # 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
66domain.set_timestepping_method('rk2')
67domain.cfl = 1.0
68domain.beta = 1.0
69#domain.limiter = "minmod_kurganov"
70#domain.limiter = "vanleer"
71#domain.limiter = "vanalbada"
72#domain.limiter = "superbee"
73#domain.limiter = "steve_minmod"
74domain.limiter = "pyvolution"
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        #xy[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        y[i] = 12.0
87    return y
88
89def elevation(x):
90    N = len(x)
91    z_infty = 10.0
92    z = zeros(N,Float)
93    L_x = 2500.0
94    A0 = 0.5*L_x
95    omega = sqrt(2*g*z_infty)/L_x
96    for i in range(N):
97        z[i] = z_infty*(x[i]**2/L_x**2)
98    return z
99
100def height(x):
101    z_infty = 10.0       ## max equilibrium water depth at lowest point.
102    L_x = 2500.0         ## width of channel
103    A0 = 0.5*L_x                  ## determines amplitudes of oscillations
104    omega = sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation
105    t=0.0
106    y = zeros(len(x),Float)
107    for i in range(len(x)):
108        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)
109    return y
110
111   
112domain.set_quantity('stage', stage)
113domain.set_quantity('elevation',elevation)
114#domain.set_quantity('height',height)
115domain.set_boundary({'exterior': Reflective_boundary(domain)})
116 
117import time
118t0 = time.time()
119#finaltime = 1122.0*0.75
120yieldstep = finaltime = 10.0
121#yieldstep = 10.0
122#plot1 = newLinePlot("Stage")
123#plot2 = newLinePlot("Xmomentum")
124C = domain.centroids
125X = domain.vertices
126StageQ = domain.quantities['stage'].vertex_values
127XmomQ = domain.quantities['xmomentum'].vertex_values
128#Bed = domain.quantities['elevation'].vertex_values
129import time
130t0 = time.time()
131for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
132    #pass
133    domain.write_time()
134    u,h,z,z_b = analytic_cannal(X.flat,t)
135    #linePlot(plot1,X,z,X,z_b,X,StageQ)
136    #linePlot(plot2,X,u*h,X,u*h,X,XmomQ)
137    HeightQ = domain.quantities['stage'].centroid_values-domain.quantities['elevation'].centroid_values
138    u,hc,z,z_b = analytic_cannal(C,t)
139    #for k in range(N):
140    #    if hc[k] < 0.0:
141    #        hc[k] = 0.0
142    error = 1.0/(N)*sum(abs(hc-HeightQ))
143    print 'Error measured at centroids', error
144    # #raw_input()         
145    # #from Gnuplot import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
146    from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
147    # #import Gnuplot 
148    X = domain.vertices
149    u,h,z,z_b = analytic_cannal(X.flat,domain.time)
150    StageQ = domain.quantities['stage'].vertex_values
151    XmomQ = domain.quantities['xmomentum'].vertex_values
152    hold(False)
153    plot1 = subplot(211)
154    plot(X,h,X,StageQ,X,z_b)
155    #plot1.set_ylim([4,11])
156    #title('?????????')
157    xlabel('Position')
158    ylabel('Stage')
159    legend(('Analytical Solution', 'Numerical Solution', 'Bed'),
160           'upper right', shadow=True)
161    plot2 = subplot(212)
162    plot(X,u*h,X,XmomQ)
163    #plot2.set_ylim([-1,25])
164    #title('?????????????????????????')
165    xlabel('Position')
166    ylabel('Xmomentum')
167    show()
Note: See TracBrowser for help on using the repository browser.