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

Last change on this file since 5827 was 5827, checked in by sudi, 15 years ago

Adding new versions of shallow_water_domain

File size: 4.2 KB
RevLine 
[5725]1import os
2from math import sqrt, sin, cos, pi, exp
[5742]3from shallow_water_domain import *
[5725]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
[5827]33N=50
[5725]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
[5743]43domain.order = 2
44domain.set_timestepping_method('rk2')
[5741]45#domain.default_time_order = 2
[5725]46domain.cfl = 1.0
[5743]47#domain.limiter = "vanleer"
[5725]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)):
[5742]63        if 800.0 <= x[i] <=1200.0:
[5738]64            y[i]=10.0
[5725]65        else:
[5742]66            y[i]=5.0
[5725]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
[5737]109domain.set_quantity('stage',stage_perturb)
[5742]110domain.set_quantity('elevation', elevation_box)
[5725]111#domain.set_quantity('xmomentum', xmom_sincos)
[5743]112#domain.order=domain.default_order
[5725]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()
[5827]125yieldstep=10.0 #30.0
126finaltime=10.0 #20.0
[5725]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   
[5728]141        from pylab import clf,plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ion
[5725]142        #print 'Test1'
[5738]143        ion()
[5731]144        hold(False)
145        #clf()
[5725]146        #print 'test 2'
147        plot1 = subplot(211)
148        #print 'test 3'
149   
[5731]150       
151        #print 'Test4'
152
[5725]153        plot(X,ElevationQ,X,HeightQ)
[5827]154        plot1.set_ylim([-1.0,11.0])
[5725]155        xlabel('Position')
156        ylabel('Stage')
[5728]157        #legend( ('Bed Elevation', 'Numerical Solution'),
158        #   'upper right', shadow=False)
159
160       
[5725]161        #legend(('Analytical Solution', 'Numerical Solution'),
162        #   'lower right', shadow=False)
163        plot2 = subplot(212)
164        plot(X,MomentumQ)
[5738]165        #plot2.set_ylim([-1,1])
[5731]166       
[5728]167        #legend( ('Numerical Solution', 'for momentum'),
168        #   'upper right', shadow=False)
[5725]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)
[5738]177        #show()
[5728]178       
[5725]179print 'That took %.2f seconds'%(time.time()-t0)
[5738]180
181raw_input("Press any key")
[5742]182
Note: See TracBrowser for help on using the repository browser.