source: trunk/anuga_work/development/pyvolution-1d/dam_steve.py @ 7924

Last change on this file since 7924 was 3401, checked in by steve, 19 years ago

Setting up analytical solution

File size: 2.6 KB
Line 
1
2import os
3from math import sqrt, pi
4from shallow_water_1d import *
5from Numeric import allclose, array, zeros, ones, Float, take, sqrt
6from config import g, epsilon
7from analytic_dam import AnalyticDam
8
9
10h0 = 5.0
11h1 = 10.0
12
13analytical_sol = AnalyticDam(h0, h1)
14
15def newLinePlot(title='Simple Plot'):
16    import Gnuplot
17    g = Gnuplot.Gnuplot(persist=1)
18    g.title(title)
19    g('set data style linespoints') 
20    g.xlabel('x')
21    g.ylabel('y')
22    return g
23
24def linePlot(g,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    g.plot(plot1,plot2)
29    #g.plot(Gnuplot.PlotItems.Data(x1.flat,y1.flat),with="linespoints")
30    #g.plot(Gnuplot.PlotItems.Data(x2.flat,y2.flat), with="lines")
31
32debug = False
33
34print "TEST 1D-SOLUTION I"
35
36L = 2000.0     # Length of channel (m)
37N = 100        # Number of compuational cells
38cell_len = L/N # Origin = 0.0
39
40points = zeros(N+1,Float)
41for i in range(N+1):
42    points[i] = i*cell_len
43
44domain = Domain(points)
45
46def stage(x):
47    y = zeros(len(x),Float)
48    for i in range(len(x)):
49        if x[i]<=1000.0:
50            y[i] = h1
51        else:
52            y[i] = h0
53    return y
54
55domain.set_quantity('stage', stage)
56
57domain.order = 2
58domain.default_order = 2
59domain.cfl = 0.8
60#domain.beta = 0.0
61print "domain.order", domain.order
62
63if (debug == True):
64    area = domain.areas
65    for i in range(len(area)):
66        if area != 20:
67            print "Cell Areas are Wrong"
68           
69    L = domain.quantities['stage'].vertex_values
70    print "Initial Stage"
71    print L
72
73domain.set_boundary({'exterior': Reflective_boundary(domain)})
74
75X = domain.vertices
76C = domain.centroids
77plot1 = newLinePlot("Stage")
78plot2 = newLinePlot("Momentum")
79
80import time
81t0 = time.time()
82yieldstep = 1.0
83finaltime = 50.0
84for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
85    domain.write_time()
86    if t > 0.0:
87        StageQ = domain.quantities['stage'].vertex_values
88        y , my = analytical_sol(X.flat,domain.time)
89        linePlot(plot1,X,StageQ,X,y)
90        MomentumQ = domain.quantities['xmomentum'].vertex_values
91        linePlot(plot2,X,MomentumQ,X,my)
92    #raw_input('press_return')
93    #pass
94
95print 'That took %.2f seconds' %(time.time()-t0)
96
97C = domain.quantities['stage'].centroid_values
98
99if (debug == True):
100    L = domain.quantities['stage'].vertex_values
101    print "Final Stage Vertex Values"
102    print L
103    print "Final Stage Centroid Values"
104    print C
105
106#f = file('test_solution_I.out', 'w')
107#for i in range(N):
108#    f.write(str(C[i]))
109#    f.write("\n")
110#f.close
111
Note: See TracBrowser for help on using the repository browser.