source:anuga_work/development/anuga_1d/dam_h_elevation.py@5743

Last change on this file since 5743 was 5743, checked in by steve, 15 years ago
File size: 4.2 KB
Line
1import os
2from math import sqrt, sin, cos, pi, exp
3from shallow_water_domain import *
4from Numeric import zeros, Float
5#from analytic_dam_sudi import AnalyticDam
6
7h0=5.0
8h1=10.0
9
10#analytical_sol=AnalyticDam(h0,h1)
11
12"""
13def newLinePlot(title='Simple Plot'):
14    import Gnuplot
15    gg=Gnuplot.Gnuplot(persist=0)
16    gg.title(title)
17    gg('set data style linespoints')
18    gg.xlabel('x')
19    gg.ylabel('y')
20    return gg
21
22def linePlot(gg, x1, y1, x2, y2):
23    import Gnuplot
24    plot1=Gnuplot.PlotItems.Data(x1.flat, y1.flat, with="linespoints")
25    plot2=Gnuplot.PlotItems.Data(x2.flat, y2.flat, with="lines 3")
26    gg.plot(plot1, plot2)
27"""
28
29
30print "TEST 1D-SOLUTION I"
31
32L=2000.0
33N=200
34
35cell_len=L/N
36
37points=zeros(N+1, Float)
38for i in range(N+1):
39    points[i]=i*cell_len
40
41domain=Domain(points)
42
43domain.order = 2
44domain.set_timestepping_method('rk2')
45#domain.default_time_order = 2
46domain.cfl = 1.0
47#domain.limiter = "vanleer"
48
49
50
51def stage(x):
52    y=zeros(len(x), Float)
53    for i in range(len(x)):
54        if x[i]<=1000.0:
55            y[i]=10.0
56        else:
57            y[i]=5.0
58    return y
59
60def stage_perturb(x):
61    y=zeros(len(x), Float)
62    for i in range(len(x)):
63        if 800.0 <= x[i] <=1200.0:
64            y[i]=10.0
65        else:
66            y[i]=5.0
67    return y
68
69def stage_sincos(x):
70    y=zeros(len(x), Float)
71    for i in range(len(x)):
72        y[i] = 5+exp(cos(pi*x[i]*0.001))+sin(pi*x[i]*0.0005)*sin(pi*x[i]*0.0005)
73    return y
74
75
76def elevation_box(x):
77    y=zeros(len(x), Float)
78    for i in range(len(x)):
79        if 500.0 <= x[i] <= 1500.0:
80            y[i] = 2.5
81        else:
82            y[i] = 0.0
83    return y
84
85def elevation_parabol(x):
86    y=zeros(len(x), Float)
87    for i in range(len(x)):
88        if 1400.0 <= x[i] <= 1600.0:
89            y[i] = 1.25*(1+cos(0.01*pi*(x[i]-500)))
90        else:
91            y[i] = 0.0
92    return y
93
94
95def elevation_sincos(x):
96    y=zeros(len(x), Float)
97    for i in range(len(x)):
98        y[i] = sin(pi*x[i]*0.0005)*sin(pi*x[i]*0.0005)
99    return y
100
101
102
103def xmom_sincos(x):
104    y=zeros(len(x), Float)
105    for i in range(len(x)):
106        y[i] = sin(cos(pi*x[i]*0.001))
107    return y
108
109domain.set_quantity('stage',stage_perturb)
110domain.set_quantity('elevation', elevation_box)
111#domain.set_quantity('xmomentum', xmom_sincos)
112#domain.order=domain.default_order
113print "domain order", domain.order
114
115domain.set_boundary({'exterior':Reflective_boundary(domain)})
116
117X=domain.vertices
118C=domain.centroids
119#plot1x=newLinePlot("Height")
120#plot2x=newLinePlot("Momentum")
121
122
123import time
124t0=time.time()
125yieldstep=45.0 #30.0
126finaltime=45.0 #20.0
127print "integral", domain.quantities['stage'].get_integral()
128for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
129    domain.write_time()
130    print "integral", domain.quantities['stage'].get_integral()
131    if t>0.0:
132        HeightQ=domain.quantities['stage'].vertex_values
133        MomentumQ=domain.quantities['xmomentum'].vertex_values
134        ElevationQ=domain.quantities['elevation'].vertex_values
135        #h, uh=analytical_sol(X.flat, domain.time)
136        #linePlot(plot1x, X, HeightQ, X, h)
137        #linePlot(plot2x, X, MomentumQ, X, uh)
138        #print "press return"
139        #pass
140
141        from pylab import clf,plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ion
142        #print 'Test1'
143        ion()
144        hold(False)
145        #clf()
146        #print 'test 2'
147        plot1 = subplot(211)
148        #print 'test 3'
149
150
151        #print 'Test4'
152
153        plot(X,ElevationQ,X,HeightQ)
154        #plot1.set_ylim([9.99,10.01])
155        xlabel('Position')
156        ylabel('Stage')
157        #legend( ('Bed Elevation', 'Numerical Solution'),
159
160
161        #legend(('Analytical Solution', 'Numerical Solution'),
163        plot2 = subplot(212)
164        plot(X,MomentumQ)
165        #plot2.set_ylim([-1,1])
166
167        #legend( ('Numerical Solution', 'for momentum'),
169
170        xlabel('Position')
171        ylabel('Xmomentum')
172
173        file = "dam_h_"
174        #file += str(number_of_cells[i])
175        file += ".eps"
176        #savefig(file)
177        #show()
178
179print 'That took %.2f seconds'%(time.time()-t0)
180
181raw_input("Press any key")
182
Note: See TracBrowser for help on using the repository browser.