source: anuga_work/development/anuga_1d/steady_flow.py @ 6453

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

Adding new versions of shallow_water_domain

File size: 4.5 KB
Line 
1import os
2#import Gnuplot
3from math import sqrt, pi
4from shallow_water_domain 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
23
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 = 400
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
136#from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ylim,xlim, rc
137#hold(False)
138#rc('text', usetex=True)
139h_error = zeros(1,Float)
140uh_error = zeros(1,Float)
141k = 0
142yieldstep = 1.0
143finaltime = 1.0
144domain.limiter = "vanleer"
145domain.order = 2
146domain.set_timestepping_method('rk2')
147domain.cfl = 1.0
148print 'THE DOMAIN LIMITER is', domain.limiter
149import time
150t0=time.time()
151for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
152    domain.write_time()
153    print "integral", domain.quantities['stage'].get_integral()
154    if t>0.0:
155        print 'Test 1'
156        from pylab import clf,plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ion
157        #ion()
158        hold(False)
159        #rc('text', usetex=True)
160        clf()
161       
162        plot1 = subplot(211)
163        print 'Test 2'
164        plot(X,w,X,StageQ,X,ElevationQ)
165        print 'Test 3'
166        xlabel('Position')
167        ylabel('Stage (m)')
168        legend(('Analytical Solution', 'Numerical Solution', 'Channel Bed'),
169               'upper right', shadow=True)
170        plot1.set_ylim([-0.1,1.0])
171        print 'Test 4'
172        plot2 = subplot(212)
173        print 'Test 5'
174        plot(X,u*h,X,XmomQ)#/(StageQ-ElevationQ))
175        print 'Test 6'
176        xlabel('x (m)')
177        ylabel(r'X-momentum (m^2/s)')
178        plot2.set_ylim([0.297,0.301])
179        h_error[k] = 1.0/(N)*sum(abs(wc-StageC))
180        uh_error[k] = 1.0/(N)*sum(abs(uc*hc-XmomC))
181        print "h_error %.10f" %(h_error[k])
182        print "uh_error %.10f"% (uh_error[k])
183        print "h_max %.10f"%max(wc-StageC)
184        print "uh max  %.10f"%max(uc*hc-XmomC)
185        #filename = "steady_flow_"
186        #filename += domain.limiter
187        #filename += str(400)
188        #filename += ".ps"
189        #savefig(filename)
190        show()
191print 'That took %.2f seconds'%(time.time()-t0)
192raw_input("Press return key!")
Note: See TracBrowser for help on using the repository browser.