source: anuga_work/development/anuga_1d/test_discontinuous_z.py @ 7815

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

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

File size: 5.0 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=4
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 = 1
19domain.set_timestepping_method('euler')
20domain.cfl = 1.0
21domain.limiter = "minmod"
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_arbitrary2)
86domain.set_boundary({'exterior':Reflective_boundary(domain)})
87X=domain.vertices.flat
88C=domain.centroids.flat
89
90yieldstep = 1.0
91finaltime = 10.0
92
93#---------------------------------------------------------
94# Start of step
95#---------------------------------------------------------
96
97domain._order_ = domain.default_order
98       
99domain.finaltime = float(finaltime)
100
101domain.yieldtime = 0.0 # Track time between 'yields'
102
103# Initialise interval of timestep sizes (for reporting only)
104domain.min_timestep = 1.0e-10
105domain.max_timestep = 1.0
106domain.number_of_steps = 0
107domain.number_of_first_order_steps = 0
108
109
110
111# Update ghosts
112domain.update_ghosts()
113
114# Initial update of vertex and edge values
115domain.distribute_to_vertices_and_edges()
116
117# Update extrema if necessary (for reporting)
118#domain.update_extrema()
119
120# Initial update boundary values
121domain.update_boundary()
122
123# Compute fluxes across each element edge
124domain.compute_fluxes()
125
126# Update timestep to fit yieldstep and finaltime
127domain.update_timestep(yieldstep, finaltime)
128
129print 'timestep = ',domain.timestep
130
131
132# Update conserved quantities
133domain.update_conserved_quantities()
134
135# Update ghosts
136domain.update_ghosts()
137
138# Update vertex and edge values
139domain.distribute_to_vertices_and_edges()
140
141# Update boundary values
142domain.update_boundary()
143
144#------------------------------------------------------------
145# End od one step
146#------------------------------------------------------------
147
148domain.write_time()
149print "integral", domain.quantities['height'].get_integral()
150StageQ=domain.quantities['stage'].vertex_values.flat
151StageUp=domain.quantities['stage'].explicit_update
152MomentumQ=domain.quantities['xmomentum'].vertex_values.flat
153MomentumUp=domain.quantities['xmomentum'].explicit_update
154
155ElevationQ=domain.quantities['elevation'].vertex_values.flat
156VelocityQ=domain.quantities['velocity'].vertex_values.flat
157
158print "momentum", MomentumQ
159
160print "StageUpdates", StageUp
161#StageUp=domain.quantities['stage'].explicit_update
162print "MomUpdates", MomentumUp
163
164hold(False)
165plot1 = subplot(311)
166#plot(X,StageQ)
167plot(X,StageQ, X,ElevationQ)
168plot1.set_ylim([-1.0,8.0])
169xlabel('Position')
170ylabel('Stage')
171legend( ('Numerical Solution', 'Bed Elevation'), 'upper right', shadow=False)
172
173plot2 = subplot(312)
174plot(X,MomentumQ)
175#plot2.set_ylim([3.998,4.002])
176legend( ('Numerical Solution', 'for momentum'), 'upper right', shadow=False)
177xlabel('Position')
178ylabel('Xmomentum')
179
180plot3 = subplot(313)
181plot(X,VelocityQ)
182#plot3.set_ylim([-5.0,30.0])
183legend( ('Numerical Solution', 'for velocity'), 'upper right', shadow=False)
184xlabel('Position')
185ylabel('Velocity')
186
187#file = "island_"
188#file += str(finaltime)
189#file += ".png"
190#savefig(file)
191show()
192
193#print 'That took %.2f seconds'%(time.time()-t0)
194#raw_input("Press return key")
Note: See TracBrowser for help on using the repository browser.