source: anuga_work/development/anuga_1d/anuga_1d_suggestion1_2_3/dry_dam_with_nonzeroslope.py @ 7576

Last change on this file since 7576 was 7576, checked in by sudi, 14 years ago
File size: 4.1 KB
Line 
1import os
2from math import sqrt, pi
3from shallow_water_domain_suggestion1 import *
4from Numeric import allclose, array, zeros, ones, Float, take, sqrt
5from config import g, epsilon
6from numpy import sin, cos, tan, arcsin, arccos, arctan
7
8def analytical_sol(C,t):
9    h1 = 10.0       # depth upstream (m)
10    L = 2000.0      # length of stream/domain (m)
11    n = len(C)      # number of cells
12    b0 = -0.005     # bottom slope
13
14    u = zeros(n,Float)
15    h = zeros(n,Float)
16    x_prime = C-L/2.0
17    x = x_prime/cos(arctan(b0))
18   
19    for i in range(n):
20        # Calculate Analytical Solution at time t > 0
21        u3 = 2.0/3.0*sqrt(g*h1)*(1.0+ b0*t*sqrt(g/h1)+x[i]/(t*sqrt(g*h1))) 
22        h3 = h1/9.0*(2.0+0.5*b0*t*sqrt(g/h1) - x[i]/(t*sqrt(g*h1)))**2
23
24        if x[i] <= -1*sqrt(g*h1)*t:
25            u[i] = 0.0
26            h[i] = h1
27        elif x[i] <= t*sqrt(g*h1)*(2.0 + 0.5*b0*t*sqrt(g/h1)):
28            u[i] = u3
29            h[i] = (0.005*x[i] + 0.025) + h3
30        else:
31            u[i] = 0.0 
32            h[i] = 0.005*x[i] + 0.025
33    return h , u*h, u
34
35#def newLinePlot(title='Simple Plot'):
36#   import Gnuplot
37#    gg = Gnuplot.Gnuplot(persist=0)
38#    gg.terminal(postscript)
39#    gg.title(title)
40#    gg('set data style linespoints')
41#    gg.xlabel('x')
42#    gg.ylabel('y')
43#    return gg
44
45#def linePlot(gg,x1,y1,x2,y2):
46#    import Gnuplot
47#    plot1 = Gnuplot.PlotItems.Data(x1.flat,y1.flat,with="linespoints")
48#    plot2 = Gnuplot.PlotItems.Data(x2.flat,y2.flat, with="lines 3")
49#    g.plot(plot1,plot2)
50
51
52
53print "TEST 1D-SOLUTION III -- DRY BED"
54
55def stage(x):
56    h1 = 10.0
57    h0 = 0.0
58    y = zeros(len(x),Float)
59    for i in range(len(x)):
60        if x[i]<=L/2.0:
61            y[i] = h1
62        else:
63            y[i] = h0
64    return y
65
66def elevation(x):
67    y=zeros(len(x), Float)
68    for i in range(len(x)):
69        y[i] = 0.005*x[i] - 5.0
70    return y
71
72
73import time
74finaltime = 10.0
75yieldstep = 10.0
76L = 2000.0     # Length of channel (m)
77number_of_cells = [810]#,200,500,1000,2000,5000,10000,20000]
78h_error = zeros(len(number_of_cells),Float)
79uh_error = zeros(len(number_of_cells),Float)
80k = 0
81for i in range(len(number_of_cells)):
82    N = int(number_of_cells[i])
83    print "Evaluating domain with %d cells" %N
84    cell_len = L/N # Origin = 0.0
85    points = zeros(N+1,Float)
86    for j in range(N+1):
87        points[j] = j*cell_len
88       
89    domain = Domain(points)
90    domain.set_quantity('stage', stage)
91    domain.set_quantity('elevation', elevation)
92   
93    domain.set_boundary({'exterior': Reflective_boundary(domain)})
94    domain.order = 2
95    domain.set_timestepping_method('rk2')
96    domain.cfl = 1.0
97    domain.limiter = "minmod"
98
99    t0 = time.time()
100    for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
101        domain.write_time()
102
103    N = float(N)
104    StageC = domain.quantities['stage'].centroid_values
105    XmomC = domain.quantities['xmomentum'].centroid_values
106    C = domain.centroids
107    h, uh, u = analytical_sol(C,domain.time)
108    h_error[k] = 1.0/(N)*sum(abs(h-StageC))
109    uh_error[k] = 1.0/(N)*sum(abs(uh-XmomC))
110    print "h_error %.10f" %(h_error[k])
111    print "uh_error %.10f"% (uh_error[k])
112    k = k+1
113    print 'That took %.2f seconds' %(time.time()-t0)
114    X = domain.vertices
115    StageQ = domain.quantities['stage'].vertex_values
116    XmomQ = domain.quantities['xmomentum'].vertex_values
117    velQ = domain.quantities['velocity'].vertex_values
118    bedQ =domain.quantities['elevation'].vertex_values
119   
120    h, uh, u = analytical_sol(X.flat,domain.time)
121    x = X.flat
122   
123    from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
124    hold(False)
125    plot1 = subplot(211)
126    plot(x,h,'b-', x,StageQ.flat,'r-', x, bedQ,'g-')
127    plot1.set_ylim([-1,11])
128    xlabel('Position')
129    ylabel('Stage')
130    legend(('Analytical Solution', 'Numerical Solution', 'Bottom elevation'),
131           'upper right', shadow=True)
132    plot2 = subplot(212)
133    plot(x,u,x,velQ.flat)
134    plot2.set_ylim([-35,35])
135   
136    xlabel('Position')
137    ylabel('Velocity')
138   
139    file = "dry_bed_"
140    file += str(number_of_cells[i])
141    file += ".eps"
142    #savefig(file)
143    show()
144   
145print "Error in height", h_error
146print "Error in xmom", uh_error
Note: See TracBrowser for help on using the repository browser.