source: anuga_work/development/analytical_solutions/oblique_shock_01.py @ 5162

Last change on this file since 5162 was 5162, checked in by steve, 17 years ago

Updated some methods for quantity. Looks like we can use old
limiting system with larger values of beta.

File size: 4.7 KB
Line 
1"""
2Example of shallow water wave equation
3consisting of an asymetrical converging channel.
4
5Copyright 2004
6Christopher Zoppou, Stephen Roberts, Ole Nielsen, Duncan Gray
7Geoscience Australia, ANU
8
9Specific methods pertaining to the 2D shallow water equation
10are imported from shallow_water
11for use with the generic finite volume framework
12
13Conserved quantities are h, uh and vh stored as elements 0, 1 and 2 in the
14numerical vector named conserved_quantities.
15
16"""
17
18######################
19# Module imports
20#
21
22#Were these used?
23#import visualise2_chris as visualise
24#import Image, ImageGrab
25
26import sys
27from os import sep
28sys.path.append('..'+sep+'pyvolution')
29
30
31from shallow_water import Domain, Constant_height
32from shallow_water import Transmissive_boundary, Reflective_boundary,\
33     Dirichlet_boundary , Transmissive_Momentum_Set_Stage_boundary
34
35from math import sqrt, cos, sin, pi
36from mesh_factory import oblique_cross
37
38
39######################
40# Domain
41#
42leny = 30.
43lenx = 40.
44f = 5
45n = 60*f
46m = 80*f
47theta = 25
48h_bc = 1.5
49p_bc = 5.0
50
51points, elements, boundary = oblique_cross(m, n, lenx, leny, theta = theta)
52domain = Domain(points, elements, boundary)
53
54
55print 'Number of Elements = ',domain.number_of_elements
56# Order of solver
57domain.default_order=2
58domain.beta_w = 0.8
59
60# Store output
61domain.store=True
62
63# Provide file name for storing output
64domain.filename="oblique_cross_%g_%g_%g_%g_%g" % (n,m,theta,h_bc,p_bc)
65print domain.filename
66
67# Output format
68domain.format="sww" #NET.CDF binary format
69                    # "dat" for ASCII
70
71
72# Visualization smoothing
73domain.smooth=False
74domain.visualise=False
75
76#######################
77#Bed-slope and friction
78class ConstantFunctionT:
79
80    def __init__(self,value):
81        self.value = value
82
83    def __call__(self,t):
84        return self.value
85
86
87shock_hb = ConstantFunctionT(h_bc)
88shock_hh = ConstantFunctionT(h_bc)
89shock_pb = ConstantFunctionT(p_bc)
90shock_pp = ConstantFunctionT(p_bc)
91
92
93def x_slope(x, y):
94    return 0*x
95
96
97
98def shock_h(x,y):
99    n = x.shape[0]
100    w = 0*x
101    for i in range(n):
102        if x[i]<0:
103            w[i] = shock_hh(0.0)
104        else:
105            w[i] = 1.0
106    return w
107
108
109def shock_p(x,y):
110    n = x.shape[0]
111    w = 0*x
112    for i in range(n):
113        if x[i]<0:
114            w[i] = shock_pp(0.0)
115        else:
116            w[i] = 0.0
117    return w
118
119
120
121
122
123domain.set_quantity('elevation', x_slope)
124domain.set_quantity('friction', 0.0)
125
126######################
127# Boundary conditions
128#
129R = Reflective_boundary(domain)
130T = Transmissive_boundary(domain)
131D = Dirichlet_boundary([shock_hb(0.0) , shock_pb(0.0), 0.0])
132Bts = Transmissive_Momentum_Set_Stage_boundary(domain, shock_hb)
133
134domain.set_boundary({'left': D, 'right': T, 'top': R, 'bottom': R})
135
136######################
137#Initial condition
138
139domain.set_quantity('stage', shock_h )
140domain.set_quantity('xmomentum',shock_p )
141
142class SetValueWhere:
143
144    def __init__(self,quantity,value=0,xrange=0.0):
145        self.quantity =  quantity
146        self.value = value
147        self.xrange = xrange
148
149    def __call__(self,x,y):
150        n = x.shape[0]
151        w = self.quantity.centroid_values
152        for i in range(n):
153            if x[i]<self.xrange:
154                w[i] = self.value
155        return w
156
157
158Stage = domain.quantities['stage']
159Xmom =  domain.quantities['xmomentum']
160Ymom =  domain.quantities['ymomentum']
161
162#for id, face in domain.boundary:
163#    print '(',id,',',face,') ',domain.boundary[(id,face)]
164
165######################
166#Evolution
167import time
168t0 = time.time()
169for t in domain.evolve(yieldstep = 0.1, finaltime = 6.0):
170    domain.write_time()
171    id = 3399
172    print Stage.get_values(location='centroids',indices=[id])
173    print Xmom.get_values(location='centroids',indices=[id])
174    print Ymom.get_values(location='centroids',indices=[id])
175    vstage =  Stage.get_values(location='centroids',indices=[id])
176    vxmom  =  Xmom.get_values(location='centroids',indices=[id])
177    id = 12719
178    print Stage.get_values(location='centroids',indices=[id])
179    print Xmom.get_values(location='centroids',indices=[id])
180    print Ymom.get_values(location='centroids',indices=[id])
181    tclean = 2.0
182#    if t > tclean-0.09 and t < tclean+0.01:
183#        print 'Cleaning up Initial profile'
184#        setstage = SetValueWhere(quantity=Stage,value=vstage[0],xrange=10)
185#        setxmom  = SetValueWhere(quantity=Xmom,value=vxmom[0],xrange=10)
186#        Stage.set_values(setstage,location='centroids')
187#        Xmom.set_values(setxmom,location='centroids')
188#        Dnew = Dirichlet_boundary([vstage[0] , vxmom[0], 0.0])
189#        domain.set_boundary({'left': Dnew, 'right': T, 'top': R, 'bottom': R})
190
191print 'That took %.2f seconds' %(time.time()-t0)
192
193#FIXME: Compute average water depth on either side of shock and compare
194#to expected values. And also Froude numbers.
195
196
197#print "saving file?"
198#im = ImageGrab.grab()
199#im.save("ccube.eps")
Note: See TracBrowser for help on using the repository browser.