# source:anuga_work/development/shallow_water_1d/subcritical_single.py@4960

Last change on this file since 4960 was 4959, checked in by steve, 16 years ago

Copied John Jakeman's and Sudi Mungkasi's 1d shallow water code

File size: 3.5 KB
Line
1import os
2from math import sqrt, pi
3from shallow_water_1d import *
4from Numeric import allclose, array, zeros, ones, Float, take, sqrt
5from config import g, epsilon
6from analytic_dam import AnalyticDam
7
8
9h0 = 5.0
10h1 = 10.0
11
12analytical_sol = AnalyticDam(h0, h1)
13
14def newLinePlot(title='Simple Plot'):
15    import Gnuplot
16    g = Gnuplot.Gnuplot(persist=1)
17    g.title(title)
18    g('set data style linespoints')
19    g.xlabel('x')
20    g.ylabel('y')
21    return g
22
23def linePlot(g,x1,y1,x2,y2):
24    import Gnuplot
25    plot1 = Gnuplot.PlotItems.Data(x1.flat,y1.flat,with="linespoints")
26    plot2 = Gnuplot.PlotItems.Data(x2.flat,y2.flat, with="lines 3")
27    g.plot(plot1,plot2)
28    #g.plot(Gnuplot.PlotItems.Data(x1.flat,y1.flat),with="linespoints")
29    #g.plot(Gnuplot.PlotItems.Data(x2.flat,y2.flat), with="lines")
30
31def stage(x):
32    y = zeros(len(x),Float)
33    for i in range(len(x)):
34        if x[i]<=1000.0:
35            y[i] = h1
36        else:
37            y[i] = h0
38    return y
39
40
41import time
42finaltime = 30.0
43yieldstep = finaltime
44
45L = 2000.0     # Length of channel (m)
46#number_of_cells = [25,50,100,200,400,800,1600,3200,6400,12800,25600]
47number_of_cells = [20]
48h_error = zeros(len(number_of_cells),Float)
49uh_error = zeros(len(number_of_cells),Float)
50k = 0
51for i in range(len(number_of_cells)):
52    N = int(number_of_cells[i])
53    print "Evaluating domain with",N,"cells"
54    cell_len = L/N # Origin = 0.0
55    points = zeros(N+1,Float)
56    for j in range(N+1):
57        points[j] = j*cell_len
58    domain = Domain(points)
59    domain.set_quantity('stage', stage)
60    domain.set_boundary({'exterior': Reflective_boundary(domain)})
61    domain.default_order = 2
62    domain.default_time_order = 2
63    print "time order", domain.default_time_order
64    domain.cfl = 1.0
65    domain.beta = 1.0
66    domain.limiter = "steve_minmod"
67    #domain.limiter = "superbee"
68    init_integral = domain.quantities['stage'].get_integral()
69    t0 = time.time()
70    for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
71        pass
72    N = float(N)
73    assert(allclose(domain.quantities['stage'].get_integral(),init_integral))
74    StageC = domain.quantities['stage'].centroid_values
75    XmomC = domain.quantities['xmomentum'].centroid_values
76    C = domain.centroids
77    h, uh = analytical_sol(C,domain.time)
78    h_error[k] = 1.0/(N)*sum(abs(h-StageC))
79    uh_error[k] = 1.0/(N)*sum(abs(uh-XmomC))
80    print "h_error %.10f" %(h_error[k])
81    print "uh_error %.10f"% (uh_error[k])
82    k = k+1
83    print 'That took %.2f seconds' %(time.time()-t0)
84    X = domain.vertices
85    StageQ = domain.quantities['stage'].vertex_values
86    XmomQ = domain.quantities['xmomentum'].vertex_values
87    h, uh = analytical_sol(X.flat,domain.time)
88    from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot#,rc
89    #rc('text', usetex=True)
90    hold(False)
91    plot1 = subplot(211)
92    plot(X,h,X,StageQ)
93    plot1.set_ylim([4,11])
94    #title('Free Surface Elevation of a Dry Dam-Break')
95    #ylabel('Stage (m)')
96    #legend(('Analytical Solution', 'Numerical Solution'),
98    #plot2 = subplot(212)
99    #plot(X,uh,X,XmomQ)
100    #plot2.set_ylim([-1,25])
101    #title('Xmomentum Profile of a Dry Dam-Break')
102    #xlabel('x (m)')
103    #ylabel(r'X-momentum (\$m^2/s\$)')
104    #legend(('Analytical Solution', 'Numerical Solution'),