source: anuga_work/development/anuga_1d/anuga_1d_suggestion1_2_3/horizontal_dam_h_2sides.py @ 7576

Last change on this file since 7576 was 7576, checked in by sudi, 14 years ago
File size: 4.8 KB
Line 
1import os
2from math import sqrt
3from shallow_water_domain_nonwellbalanced3 import *
4#from channel_domain import *
5from Numeric import zeros, Float
6from analytic_dam_sudi import AnalyticDam
7
8h0=5.0
9h1=10.0
10
11analytical_sol=AnalyticDam(h0,h1)
12
13"""
14def newLinePlot(title='Simple Plot'):
15    import Gnuplot
16    gg=Gnuplot.Gnuplot(persist=0)
17    gg.title(title)
18    gg('set data style linespoints')
19    gg.xlabel('x')
20    gg.ylabel('y')
21    return gg
22
23def linePlot(gg, x1, y1, x2, y2):
24    import Gnuplot
25    plot1=Gnuplot.PlotItems.Data(x1.flat, y1.flat, with="linespoints")
26    plot2=Gnuplot.PlotItems.Data(x2.flat, y2.flat, with="lines 3")
27    gg.plot(plot1, plot2)
28"""
29
30
31print "TEST 1D-SOLUTION I"
32
33L=2000.0
34N=400
35
36cell_len=L/N
37
38points=zeros(N+1, Float)
39for i in range(N+1):
40    points[i]=i*cell_len
41
42domain=Domain(points)
43
44domain.order = 2
45domain.set_timestepping_method('rk2')
46domain.cfl = 1.0
47domain.limiter = "vanleer"
48
49
50
51def height(x):
52    y=zeros(len(x), Float)
53    for i in range (len(x)):
54        if x[i]<=L/4.0:
55            y[i]=0.0 #h0
56        elif x[i]<=3*L/4.0:
57            y[i]=h1
58        else:
59            y[i]=h0
60    return y
61
62domain.set_quantity('stage',height)
63print "domain order=", domain.order
64print "domain limiter=", domain.limiter
65
66domain.set_boundary({'exterior':Reflective_boundary(domain)})
67
68X=domain.vertices
69C=domain.centroids
70#plot1x=newLinePlot("Height")
71#plot2x=newLinePlot("Momentum")
72
73import time
74yieldstep=finaltime=10.0
75t0=time.time()
76
77"""
78for i in range(10):
79    L=2000.0
80    N=400
81    cell_len=L/N
82    points=zeros(N+1, Float)
83    for i in range(N+1):
84        points[i]=i*cell_len
85    domain=Domain(points)
86    domain.order = 2
87    domain.set_timestepping_method('rk2')
88    domain.cfl = 1.0
89    domain.limiter = "vanleer"
90    domain.set_quantity('stage',height)
91    domain.set_boundary({'exterior':Reflective_boundary(domain)})
92    #print "integral", domain.quantities['stage'].get_integral()
93    for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
94        domain.write_time()
95    #print "integral", domain.quantities['stage'].get_integral()
96    #print '==================================================================='
97print 'The average time is %.2f seconds'%((time.time()-t0)/10.0)
98"""
99while finaltime<10.1:
100    yieldstep=finaltime
101    for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
102        domain.write_time()
103    if t>0.0:
104        N = float(N)
105        StageC = domain.quantities['stage'].centroid_values
106        XmomC = domain.quantities['xmomentum'].centroid_values
107        VelC = domain.quantities['velocity'].centroid_values
108        C = domain.centroids
109        hC, uhC, uC = analytical_sol(C,domain.time)
110        h_error = cell_len*sum(abs(hC-StageC))
111        uh_error = cell_len*sum(abs(uhC-XmomC))
112        u_error = cell_len*sum(abs(uC-VelC))
113        print "h_error %.10f" %(h_error)
114        print "uh_error %.10f"%(uh_error)
115        print "u_error %.10f" %(u_error)
116        print 'That took %.2f seconds' %(time.time()-t0)
117        X = domain.vertices
118        StageQ = domain.quantities['stage'].vertex_values
119        XmomQ = domain.quantities['xmomentum'].vertex_values
120        VelQ = domain.quantities['velocity'].vertex_values
121
122        h, uh, u = analytical_sol(X.flat, domain.time)
123        #linePlot(plot1x, X, HeightQ, X, h)
124        #linePlot(plot2x, X, MomentumQ, X, uh)
125        #print "press return"
126        #pass
127   
128        from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
129        hold(False)
130       
131        plot1 = subplot(311)
132        plot(X,h,X,StageQ)
133        plot1.set_ylim([-1,11])
134        plot1.set_xlim([-100.0,2000.0])
135        #xlabel('Position')
136        ylabel('Stage')
137        #legend(('Analytical Solution', 'Numerical Solution'),
138        #   'lower right', shadow=False)
139       
140        plot2 = subplot(312)
141        plot(X,uh,X,XmomQ)
142        #plot2.set_ylim([-10,1])
143        plot2.set_xlim([-100.0,2000.0])
144        #legend(('Analytical Solution', 'Numerical Solution'),
145        #        'lower right', shadow=False)
146        #xlabel('Position')
147        ylabel('Momentum')
148
149       
150        plot3 = subplot(313)
151        plot(X,u,X,VelQ)
152        #plot3.set_ylim([-5,35])
153        plot3.set_xlim([-100.0,2000.0])
154        legend(('Analytical Solution', 'Numerical Solution'),
155                'lower right', shadow=False)
156        xlabel('Position')
157        ylabel('Velocity')
158       
159       
160        #file = "dam_h_"
161        #file += str(number_of_cells[i])
162        #file += ".eps"
163        #savefig(file)
164        #show()
165
166        #file = "dam_"
167        #file += str(finaltime)
168        #file += ".png"
169        #savefig(file)
170        #show()
171    print "The domain.limiter is", domain.limiter
172    print 'That took %.2f seconds'%(time.time()-t0)
173    print '============================================================================='
174    finaltime = finaltime + 5.0
Note: See TracBrowser for help on using the repository browser.