# source:anuga_work/development/anuga_1d/discontinuous_z.py@6713

Last change on this file since 6713 was 6713, checked in by steve, 14 years ago

Fixed pylab so that it would work on bogong, by flattening the arrays for plotting

File size: 3.9 KB
Line
1import os
2import time
3from math import sqrt, sin, cos, pi, exp
4from shallow_water_domain_suggestion2 import *
5from Numeric import zeros, Float
6from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
7
8print "TEST 1D-SOLUTION I"
9L=2000.0
10N=200
11
12cell_len=L/N
13points=zeros(N+1, Float)
14for i in range(N+1):
15    points[i]=i*cell_len
16
17domain=Domain(points)
18domain.order = 2
19domain.set_timestepping_method('rk2')
20domain.cfl = 1.0
21domain.limiter = "vanleer"
22
23def stage_flat(x):
24    y=zeros(len(x), Float)
25    for i in range(len(x)):
26        y[i]=4.0
27    return y
28
29def elevation_arbitrary(x):
30    y=zeros(len(x), Float)
31    for i in range(len(x)):
32        if 0 <= x[i] < 200.0:
33            y[i] = -0.01*(x[i]-200) + 4.0
34        elif 200.0 <= x[i] < 300.0:
35            y[i] = -0.02*(x[i]-200) + 4.0
36        elif 300.0 <= x[i] < 400.0:
37            y[i] = -0.01*(x[i]-300) + 2.0
38        elif 400.0 <= x[i] < 550.0:
39            y[i] = (-1/75.0)*(x[i]-400.0) + 2.0
40        elif 550.0 <= x[i] < 700.0:
41            y[i] = (1/11250)*(x[i]-550)*(x[i]-550)
42        elif 700.0 <= x[i] < 800.0:
43            y[i] = 0.03*(x[i]-700)
44        elif 800.0 <= x[i] < 900.0:
45            y[i] = -0.03*(x[i]-800) + 3.0
46        elif 900.0 <= x[i] < 1000.0:
47            y[i] = 6.0
48        elif 1000.0 <= x[i] < 1400.0:
49            y[i] = (-1.0/20000)*(x[i]-1000)*(x[i]-1400)
50        elif 1400.0 <= x[i] < 1500.0:
51            y[i] = 0.0
52        elif 1500.0 <= x[i] < 1700.0:
53            y[i] = 3.0
54        elif 1700.0 <= x[i] < 1800.0:
55            y[i] = -0.03*(x[i]-1700) + 3.0
56        else:
57            y[i] = (3.0/40000)*(x[i]-1800)*(x[i]-1800) + 2.0
58    return y
59
60def elevation_arbitrary1(x):
61    y=zeros(len(x), Float)
62    for i in range(len(x)):
63        if 0 <= x[i] < 200.0:
64            y[i] = -0.01*(x[i]-200) + 5.0
65        elif 200.0 <= x[i] < 900.0:
66            y[i] = 0.0
67        elif 900.0 <= x[i] < 1000.0:    # This is the island
68            y[i] = 5.0
69        elif 1000.0 <= x[i] < 1800.0:
70            y[i] = 0.0
71        else:
72            y[i] = 0.03*(x[i]-1700) + 3.0
73    return y
74
75def elevation_arbitrary2(x):
76    y=zeros(len(x), Float)
77    for i in range(len(x)):
78        if 0 <= x[i] < 1000.0:          # This is the island
79            y[i] = 5.0
80        else:
81            y[i] = 0.0
82    return y
83
84domain.set_quantity('stage',stage_flat)
85domain.set_quantity('elevation',elevation_arbitrary)
86domain.set_boundary({'exterior':Reflective_boundary(domain)})
87X=domain.vertices.flat
88C=domain.centroids.flat
89
90t0=time.time()
91yieldstep=finaltime=60.0
92while finaltime < 60.2:
93    for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
94        domain.write_time()
95        print "t=",t
96        print "integral", domain.quantities['height'].get_integral()
97        StageQ=domain.quantities['stage'].vertex_values.flat
98        MomentumQ=domain.quantities['xmomentum'].vertex_values.flat
99        ElevationQ=domain.quantities['elevation'].vertex_values.flat
100        VelocityQ=domain.quantities['velocity'].vertex_values.flat
101
102        hold(False)
103        plot1 = subplot(311)
104        plot(X,StageQ, X,ElevationQ)
105        plot1.set_ylim([-1.0,8.0])
106        xlabel('Position')
107        ylabel('Stage')
108        legend( ('Numerical Solution', 'Bed Elevation'), 'upper right', shadow=False)
109
110        plot2 = subplot(312)
111        plot(X,MomentumQ)
112        #plot2.set_ylim([3.998,4.002])
113        legend( ('Numerical Solution', 'for momentum'), 'upper right', shadow=False)
114        xlabel('Position')
115        ylabel('Xmomentum')
116
117        plot3 = subplot(313)
118        plot(X,VelocityQ)
119        #plot2.set_ylim([-5.0,30.0])
120        legend( ('Numerical Solution', 'for velocity'), 'upper right', shadow=False)
121        xlabel('Position')
122        ylabel('Velocity')
123
124        file = "island_"
125        file += str(finaltime)
126        file += ".png"
127        #savefig(file)
128        show()
129
130    print 'That took %.2f seconds'%(time.time()-t0)
131    finaltime = finaltime + 10.0
132raw_input("Press return key")
Note: See TracBrowser for help on using the repository browser.