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

Last change on this file since 6454 was 5587, checked in by steve, 16 years ago

File size: 2.9 KB
Line
1import os
2from math import sqrt
3#from shallow_water_h import *
4from shallow_water_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.default_order = 2
45domain.default_time_order = 1
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) #('height', height)
63domain.order=domain.default_order
64print "domain order", domain.order
65
66domain.set_boundary({'exterior':Reflective_boundary(domain)})
67
68X=domain.vertices
69C=domain.centroids
70#plot1x=newLinePlot("Height")
71#plot2x=newLinePlot("Momentum")
72
73
74import time
75t0=time.time()
76yieldstep=30.0
77finaltime=20.0
78print "integral", domain.quantities['stage'].get_integral() #['height'].get_integral()
79for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
80    domain.write_time()
81    print "integral", domain.quantities['stage'].get_integral() #['height'].get_integral()
82    if t>0.0:
83        HeightQ=domain.quantities['stage'].vertex_values #['height'].vertex_values
84        MomentumQ=domain.quantities['xmomentum'].vertex_values
85        h, uh=analytical_sol(X.flat, domain.time)
86        #linePlot(plot1x, X, HeightQ, X, h)
87        #linePlot(plot2x, X, MomentumQ, X, uh)
88        #print "press return"
89        #pass
90
91        from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
92        #print 'Test1'
93        hold(False)
94        #print 'test 2'
95        plot1 = subplot(211)
96        #print 'test 3'
97
98        plot(X,h,X,HeightQ)
99        #print 'Test4'
100        plot1.set_ylim([0,11])
101        xlabel('Position')
102        ylabel('Stage')
103        #legend(('Analytical Solution', 'Numerical Solution'),
105        plot2 = subplot(212)
106        plot(X,uh,X,MomentumQ)
107        #plot2.set_ylim([-5,35])
108        legend(('Analytical Solution', 'Numerical Solution'),