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

Last change on this file since 5565 was 5565, checked in by steve, 14 years ago

Added profiling script (which shows a speedup of 3 using C code for compute fluxes)

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