source:anuga_work/development/shallow_water_1d/dry_dam.py@4960

Last change on this file since 4960 was 4960, checked in by steve, 16 years ago
File size: 4.3 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
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-L/2.0
19    #x = zeros(n,Float)
20    #for i in range(n):
21    #    x[i] = C[i]-1000.0
22
23    for i in range(n):
24        # Calculate Analytical Solution at time t > 0
25        u3 = 2.0/3.0*(sqrt(g*h1)+x[i]/t)
26        h3 = 4.0/(9.0*g)*(sqrt(g*h1)-x[i]/(2.0*t))*(sqrt(g*h1)-x[i]/(2.0*t))
27
28        if ( x[i] <= -t*sqrt(g*h1) ):
29            u[i] = 0.0
30            h[i] = h1
31        elif ( x[i] <= 2.0*t*sqrt(g*h1) ):
32            u[i] = u3
33            h[i] = h3
34        else:
35            u[i] = 0.0
36            h[i] = h0
37
38    return h , u*h
39
40def newLinePlot(title='Simple Plot'):
41    import Gnuplot
42    gg = Gnuplot.Gnuplot(persist=0)
43    gg.terminal(postscript)
44    gg.title(title)
45    gg('set data style linespoints')
46    gg.xlabel('x')
47    gg.ylabel('y')
48    return gg
49
50def linePlot(gg,x1,y1,x2,y2):
51    import Gnuplot
52    plot1 = Gnuplot.PlotItems.Data(x1.flat,y1.flat,with="linespoints")
53    plot2 = Gnuplot.PlotItems.Data(x2.flat,y2.flat, with="lines 3")
54    g.plot(plot1,plot2)
55
56#debug = False
57
58print "TEST 1D-SOLUTION III -- DRY BED"
59
60L = 2000.0     # Length of channel (m)
61N = 100        # Number of compuational cells
62cell_len = L/N # Origin = 0.0
63
64points = zeros(N+1,Float)
65for i in range(N+1):
66    points[i] = i*cell_len
67
68domain = Domain(points)
69
70def stage(x):
71    y = zeros(len(x),Float)
72    for i in range(len(x)):
73        if x[i]<=1000.0:
74            y[i] = 10.0
75        else:
76            y[i] = 0.0
77    return y
78
79
80import time
81
82finaltime = 30.0
83yieldstep = 1
84L = 2000.0     # Length of channel (m)
85number_of_cells = [100]#,200,500,1000,2000,5000,10000,20000]
86h_error = zeros(len(number_of_cells),Float)
87uh_error = zeros(len(number_of_cells),Float)
88k = 0
89for i in range(len(number_of_cells)):
90    N = int(number_of_cells[i])
91    print "Evaluating domain with %d cells" %N
92    cell_len = L/N # Origin = 0.0
93    points = zeros(N+1,Float)
94    for j in range(N+1):
95        points[j] = j*cell_len
96
97    domain = Domain(points)
98
99    domain.set_quantity('stage', stage)
100    domain.set_boundary({'exterior': Reflective_boundary(domain)})
101    domain.default_order = 2
102    domain.default_time_order = 1
103    domain.cfl = 1.0
104    domain.limiter = "vanleer"
105
106    t0 = time.time()
107    for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
108        pass
109    N = float(N)
110    StageC = domain.quantities['stage'].centroid_values
111    XmomC = domain.quantities['xmomentum'].centroid_values
112    C = domain.centroids
113    h, uh = analytical_sol(C,domain.time)
114    h_error[k] = 1.0/(N)*sum(abs(h-StageC))
115    uh_error[k] = 1.0/(N)*sum(abs(uh-XmomC))
116    print "h_error %.10f" %(h_error[k])
117    print "uh_error %.10f"% (uh_error[k])
118    k = k+1
119    print 'That took %.2f seconds' %(time.time()-t0)
120    X = domain.vertices.flat
121    StageQ = domain.quantities['stage'].vertex_values.flat
122    XmomQ = domain.quantities['xmomentum'].vertex_values.flat
123    h, uh = analytical_sol(X,domain.time)
124    #from pylab import * # this command stuffs up setting of new domain
125    from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
126    print 'test2'
127    #hold(False)
128    #subplot(211)
129    #plot1 = plot(X,h,X,StageQ)
130    hold(False)
131    plot1 = subplot(211)
132    plot(X,h,X,StageQ)
133    plot1.set_ylim([-1,11])
134    #plot1 = plot(X,StageQ,'o',X,h)
135    #set(markersize=1)
136    #title('Free Surface Elevation of a Dry Dam-Break')
137    xlabel('Position')
138    ylabel('Stage')
139    legend(('Analytical Solution', 'Numerical Solution'),
141    plot2 = subplot(212)
142    plot(X,uh,X,XmomQ)
143    plot2.set_ylim([-5,35])
144    #title('Xmomentum Profile of a Dry Dam-Break')
145    xlabel('Position')
146    ylabel('Xmomentum')
147    #legend(('Analytical Solution', 'Numerical Solution'),