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: 3.9 KB
Line
1import os
2import Gnuplot
3from math import sqrt, pi
4from shallow_water_1d import *
5from Numeric import allclose, array, zeros, ones, Float, take, sqrt
6from 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
19h_0 = 0.5
20u_0 = 0.6
21h_1 = 0.5
22u_1 = 0.6
23g = 9.81
24
25#h_c = (q**2/g)**(1/3)
26#u_c = sqrt(g*h_c)
27
28q = 0.3
29g = 9.81
30
31b = 4.0
32xmax = 25.0
33z_bmax = 0.2
34Fr_0 = u_0/sqrt(g*h_0)
35
36def stage(x):
37    y,u,h = analytical_sol(x)
38    return y
39
40
41def 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
54def xmomentum(x):
55    return u_0*h_0
56
57
58def 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
73def f(D,H):
74    return D**3+D**2*(H-1.0-Fr_0**2/2.0)+Fr_0**2/2.0
75
76def 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
80def 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
99Fr_0 = u_0/sqrt(g*h_0)
100zb=0.2
101H = zb/h_0
102D = bisection(f,1,0.2/0.5,H)
103
104L = 50.0
105N = 50 #800
106cell_len = L/N
107points = zeros(N+1,Float)
108for j in range(N+1):
109    points[j] = j*cell_len
110boundary = { (0,0): 'left',(N-1,1): 'right'}
111domain = Domain(points,boundary)
112D1 = Dirichlet_boundary([h_0,u_0*h_0])
113D2 = Dirichlet_boundary([h_1,u_1*h_1])
114domain.set_boundary({'left':D1,'right':D2})
115domain.set_quantity('elevation',elevation)
116domain.set_quantity('stage',stage)
117domain.set_quantity('xmomentum',xmomentum)
118X = domain.vertices
119C = domain.centroids
120
121Stage = domain.quantities['stage']
122Xmom = domain.quantities['xmomentum']
123
124StageQ = Stage.vertex_values
125StageC = Stage.centroid_values
126XmomQ = Xmom.vertex_values
127XmomC = Xmom.centroid_values
128ElevationQ = domain.quantities['elevation'].vertex_values
129
130stage_bdry = Stage.boundary_values
131xmom_bdry =  Xmom.boundary_values
132
133w,u,h = analytical_sol(X.flat)
134wc, uc, hc = analytical_sol(C)
135
136from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ylim,xlim, rc
137hold(False)
138rc('text', usetex=True)
139h_error = zeros(1,Float)
140uh_error = zeros(1,Float)
141k = 0
142yieldstep = 25.0
143finaltime = 25.0
144domain.limiter = "vanleer"
145#domain.limiter = "steve_minmod"
146domain.default_order = 2
147domain.default_time_order = 2#check order is not working
148for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
149    domain.write_time()
150print 'Test 1'
151plot1 = subplot(211)
152print 'Test 2'
153plot(X,w,X,StageQ,X,ElevationQ)
154print 'Test 3'
155#xlabel('Position')
156#ylabel('Stage (m)')
157#legend(('Analytical Solution', 'Numerical Solution', 'Channel Bed'),
159plot1.set_ylim([-0.1,1.0])
160print 'Test 4'
161plot2 = subplot(212)
162print 'Test 5'
163plot(X,u*h,X,XmomQ)#/(StageQ-ElevationQ))
164print 'Test 6'
165#xlabel('x (m)')
166#ylabel(r'X-momentum (\$m^2/s\$)')
167h_error[k] = 1.0/(N)*sum(abs(wc-StageC))
168uh_error[k] = 1.0/(N)*sum(abs(uc*hc-XmomC))
169print "h_error %.10f" %(h_error[k])
170print "uh_error %.10f"% (uh_error[k])
171print "h_max %.10f"%max(wc-StageC)
172print "uh max  %.10f"%max(uc*hc-XmomC)