source: trunk/anuga_work/development/shallow_water_1d_old/dam_h_sudi.py @ 7884

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

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

File size: 2.8 KB
Line 
1import os
2from math import sqrt
3from shallow_water_h import *
4from Numeric import zeros, Float
5from analytic_dam_sudi import AnalyticDam
6
7h0=5.0
8h1=10.0
9
10analytical_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=100
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 = 1
44domain.default_time_order = 1
45#domain.cfl = 1.0
46#domain.limiter = "minmod"
47
48
49
50
51
52   
53   
54def height(x):
55    y=zeros(len(x), Float)
56    for i in range (len(x)):
57        if x[i]<=L/4.0:
58            y[i]=0.0 #h0
59        elif x[i]<=3*L/4.0:
60            y[i]=h1
61        else:
62            y[i]=h0
63    return y
64
65domain.set_quantity('height', height)
66domain.order=domain.default_order
67print "domain order", domain.order
68
69domain.set_boundary({'exterior':Reflective_boundary(domain)})
70
71X=domain.vertices
72C=domain.centroids
73#plot1x=newLinePlot("Height")
74#plot2x=newLinePlot("Momentum")
75
76
77import time
78t0=time.time()
79yieldstep=30.0
80finaltime=30.0
81print "integral", domain.quantities['height'].get_integral()
82for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
83    domain.write_time()
84    print "integral", domain.quantities['height'].get_integral()
85    if t>0.0:
86        HeightQ=domain.quantities['height'].vertex_values
87        MomentumQ=domain.quantities['xmomentum'].vertex_values
88        h, uh=analytical_sol(X.flat, domain.time)
89        #linePlot(plot1x, X, HeightQ, X, h)
90        #linePlot(plot2x, X, MomentumQ, X, uh)
91        #print "press return"
92        #pass
93   
94        from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
95        #print 'Test1'
96        hold(False)
97        #print 'test 2'
98        plot1 = subplot(211)
99        #print 'test 3'
100   
101        plot(X,h,X,HeightQ)
102        #print 'Test4'
103        plot1.set_ylim([0,11])
104        xlabel('Position')
105        ylabel('Stage')
106        #legend(('Analytical Solution', 'Numerical Solution'),
107        #   'lower right', shadow=False)
108        plot2 = subplot(212)
109        plot(X,uh,X,MomentumQ)
110        #plot2.set_ylim([-5,35])
111        legend(('Analytical Solution', 'Numerical Solution'),
112           'lower right', shadow=False)
113   
114        xlabel('Position')
115        ylabel('Xmomentum')
116       
117        file = "dam_h_"
118        #file += str(number_of_cells[i])
119        file += ".eps"
120        #savefig(file)
121        show()
122
123print 'That took %.2f seconds'%(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.