source: anuga_work/development/anuga_1d/dam_h_sudi.py @ 7818

Last change on this file since 7818 was 6695, checked in by sudi, 14 years ago

Just modifying.

File size: 2.8 KB
Line 
1import os
2from math import sqrt
3#from shallow_water_h import *
4from shallow_water_domain_suggestion2 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=200
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)
63
64domain.set_boundary({'exterior':Reflective_boundary(domain)})
65
66X=domain.vertices
67C=domain.centroids
68#plot1x=newLinePlot("Height")
69#plot2x=newLinePlot("Momentum")
70
71
72import time
73t0=time.time()
74yieldstep=20.0
75finaltime=20.0
76print "integral", domain.quantities['stage'].get_integral() #['height'].get_integral()
77for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
78    domain.write_time()
79    print "integral", domain.quantities['stage'].get_integral() #['height'].get_integral()
80    if t>0.0:
81        HeightQ=domain.quantities['stage'].vertex_values #['height'].vertex_values
82        MomentumQ=domain.quantities['xmomentum'].vertex_values
83        h, uh=analytical_sol(X.flat, domain.time)
84        #linePlot(plot1x, X, HeightQ, X, h)
85        #linePlot(plot2x, X, MomentumQ, X, uh)
86        #print "press return"
87        #pass
88   
89        from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
90        #print 'Test1'
91        hold(False)
92        #print 'test 2'
93        plot1 = subplot(211)
94        #print 'test 3'
95   
96        plot(X,h,X,HeightQ)
97        #print 'Test4'
98        plot1.set_ylim([0,11])
99        xlabel('Position')
100        ylabel('Stage')
101        #legend(('Analytical Solution', 'Numerical Solution'),
102        #   'lower right', shadow=False)
103        plot2 = subplot(212)
104        plot(X,uh,X,MomentumQ)
105        #plot2.set_ylim([-5,35])
106        legend(('Analytical Solution', 'Numerical Solution'),
107           'lower right', shadow=False)
108   
109        xlabel('Position')
110        ylabel('Xmomentum')
111       
112        file = "dam_h_"
113        #file += str(number_of_cells[i])
114        file += ".eps"
115        #savefig(file)
116        show()
117
118print 'That took %.2f seconds'%(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.