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

Last change on this file since 5725 was 5725, checked in by sudi, 16 years ago
File size: 4.0 KB
Line 
1import os
2from math import sqrt, sin, cos, pi, exp
3from shallow_water_domain_new 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=800
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.default_order = 2
44domain.default_time_order = 2
45domain.cfl = 1.0
46domain.limiter = "vanleer"
47
48
49
50def stage(x):
51    y=zeros(len(x), Float)
52    for i in range(len(x)):
53        if x[i]<=1000.0:
54            y[i]=10.0
55        else:
56            y[i]=5.0
57    return y
58
59def stage_perturb(x):
60    y=zeros(len(x), Float)
61    for i in range(len(x)):
62        if 1100.0 <= x[i] <=1200.0:
63            y[i]=10.0
64        else:
65            y[i]=10.0
66    return y
67
68def stage_sincos(x):
69    y=zeros(len(x), Float)
70    for i in range(len(x)):
71        y[i] = 5+exp(cos(pi*x[i]*0.001))+sin(pi*x[i]*0.0005)*sin(pi*x[i]*0.0005)
72    return y
73
74
75def elevation_box(x):
76    y=zeros(len(x), Float)
77    for i in range(len(x)):
78        if 500.0 <= x[i] <= 1500.0:
79            y[i] = 2.5
80        else:
81            y[i] = 0.0
82    return y
83
84def elevation_parabol(x):
85    y=zeros(len(x), Float)
86    for i in range(len(x)):
87        if 1400.0 <= x[i] <= 1600.0:
88            y[i] = 1.25*(1+cos(0.01*pi*(x[i]-500)))
89        else:
90            y[i] = 0.0
91    return y
92
93
94def elevation_sincos(x):
95    y=zeros(len(x), Float)
96    for i in range(len(x)):
97        y[i] = sin(pi*x[i]*0.0005)*sin(pi*x[i]*0.0005)
98    return y
99
100
101
102def xmom_sincos(x):
103    y=zeros(len(x), Float)
104    for i in range(len(x)):
105        y[i] = sin(cos(pi*x[i]*0.001))
106    return y
107
108domain.set_quantity('stage',stage_perturb)
109domain.set_quantity('elevation', elevation_parabol)
110#domain.set_quantity('xmomentum', xmom_sincos)
111domain.order=domain.default_order
112print "domain order", domain.order
113
114domain.set_boundary({'exterior':Reflective_boundary(domain)})
115
116X=domain.vertices
117C=domain.centroids
118#plot1x=newLinePlot("Height")
119#plot2x=newLinePlot("Momentum")
120
121
122import time
123t0=time.time()
124yieldstep=0.5 #30.0
125finaltime=10.0 #20.0
126print "integral", domain.quantities['stage'].get_integral()
127for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
128    domain.write_time()
129    print "integral", domain.quantities['stage'].get_integral()
130    if t>0.0:
131        HeightQ=domain.quantities['stage'].vertex_values
132        MomentumQ=domain.quantities['xmomentum'].vertex_values
133        ElevationQ=domain.quantities['elevation'].vertex_values
134        #h, uh=analytical_sol(X.flat, domain.time)
135        #linePlot(plot1x, X, HeightQ, X, h)
136        #linePlot(plot2x, X, MomentumQ, X, uh)
137        #print "press return"
138        #pass
139   
140        from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ion
141        #print 'Test1'
142        hold(False)
143        #print 'test 2'
144        plot1 = subplot(211)
145        #print 'test 3'
146   
147        plot(X,ElevationQ,X,HeightQ)
148        #print 'Test4'
149        plot1.set_ylim([-1,12])
150        xlabel('Position')
151        ylabel('Stage')
152        legend( ('Bed Elevation', 'Numerical Solution'),
153           'upper right', shadow=False)
154        #legend(('Analytical Solution', 'Numerical Solution'),
155        #   'lower right', shadow=False)
156        plot2 = subplot(212)
157        plot(X,MomentumQ)
158        #plot2.set_ylim([-5,35])
159        legend( ('Numerical Solution', 'for momentum'),
160           'upper right', shadow=False)
161   
162        xlabel('Position')
163        ylabel('Xmomentum')
164       
165        file = "dam_h_"
166        #file += str(number_of_cells[i])
167        file += ".eps"
168        #savefig(file)
169        show()
170
171print 'That took %.2f seconds'%(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.