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

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