source: anuga_work/development/anuga_1d/dry_dam_sudi.py @ 5832

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