source: anuga_work/development/shallow_water_1d/dam_h.py @ 4959

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

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

File size: 2.1 KB
Line 
1import os
2#from math import sqrt, pi
3from shallow_water_h import Domain, Reflective_boundary
4from Numeric import zeros, Float
5#from config import g, epsilon
6from analytic_dam import AnalyticDam
7
8
9h0 = 5.0
10h1 = 10.0
11
12analytical_sol = AnalyticDam(h0, h1)
13
14def newLinePlot(title='Simple Plot'):
15    import Gnuplot
16    g = Gnuplot.Gnuplot(persist=0)
17    g.title(title)
18    g('set data style linespoints') 
19    g.xlabel('x')
20    g.ylabel('y')
21    return g
22
23def linePlot(g,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    g.plot(plot1,plot2)
28   
29
30print "TEST 1D-SOLUTION I"
31
32L = 2000.0     # Length of channel (m)
33N = 100        # Number of computational cells
34cell_len = L/N # Origin = 0.0
35
36points = zeros(N+1,Float)
37for i in range(N+1):
38    points[i] = i*cell_len
39
40domain = Domain(points)
41
42def height(x):
43    y = zeros(len(x),Float)
44    for i in range(len(x)):
45        if x[i]<=1000.0:
46            y[i] = h1
47        else:
48            y[i] = h0
49    return y
50
51
52domain.set_quantity('height', height)
53
54domain.default_order = 2
55domain.order = 2
56domain.default_time_order = 1
57domain.cfl = 0.8
58#domain.beta = 1.0
59domain.split = False
60print "domain.order", domain.order
61domain.limiter = "minmod"
62
63
64domain.set_boundary({'exterior': Reflective_boundary(domain)})
65
66X = domain.vertices
67C = domain.centroids
68plot1 = newLinePlot("Height")
69plot2 = newLinePlot("Momentum")
70
71import time
72t0 = time.time()
73yieldstep = 1.0
74finaltime = 20.0
75print "integral", domain.quantities['height'].get_integral()
76for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
77    domain.write_time()
78    print "integral", domain.quantities['height'].get_integral()
79    if t > 0.0:
80        HeightQ = domain.quantities['height'].vertex_values
81        y , my = analytical_sol(X.flat,domain.time)
82        linePlot(plot1,X,HeightQ,X,y)
83        MomentumQ = domain.quantities['xmomentum'].vertex_values
84        linePlot(plot2,X,MomentumQ,X,my)
85    #raw_input('press_return')
86
87    #pass
88
89print 'That took %.2f seconds' %(time.time()-t0)
90
91#C = domain.quantities['height'].centroid_values
92
Note: See TracBrowser for help on using the repository browser.