source: anuga_work/development/anuga_1d/anuga_1d_suggestion1_2_3/test_discontinuous_z.py @ 7576

Last change on this file since 7576 was 7576, checked in by sudi, 14 years ago
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
88C=domain.centroids
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
151StageUp=domain.quantities['stage'].explicit_update
152MomentumQ=domain.quantities['xmomentum'].vertex_values
153MomentumUp=domain.quantities['xmomentum'].explicit_update
154
155ElevationQ=domain.quantities['elevation'].vertex_values
156VelocityQ=domain.quantities['velocity'].vertex_values
157
158print "momentum", MomentumQ
159
160print "StageUpdates", StageUp
161#StageUp=domain.quantities['stage'].explicit_update
162print "MomUpdates", MomentumUp
163hold(False)
164plot1 = subplot(311)
165plot(X,StageQ, X,ElevationQ)
166plot1.set_ylim([-1.0,8.0])
167xlabel('Position')
168ylabel('Stage')
169legend( ('Numerical Solution', 'Bed Elevation'), 'upper right', shadow=False)
170
171plot2 = subplot(312)
172plot(X,MomentumQ)
173#plot2.set_ylim([3.998,4.002])
174legend( ('Numerical Solution', 'for momentum'), 'upper right', shadow=False)
175xlabel('Position')
176ylabel('Xmomentum')
177
178plot3 = subplot(313)
179plot(X,VelocityQ)
180#plot2.set_ylim([-5.0,30.0])
181legend( ('Numerical Solution', 'for velocity'), 'upper right', shadow=False)
182xlabel('Position')
183ylabel('Velocity')
184
185#file = "island_"
186#file += str(finaltime)
187#file += ".png"
188#savefig(file)
189show()
190
191#print 'That took %.2f seconds'%(time.time()-t0)
192#raw_input("Press return key")
Note: See TracBrowser for help on using the repository browser.