source: inundation/analytical solutions/Analytical solution_oblique_shock_01.py @ 2152

Last change on this file since 2152 was 1949, checked in by steve, 19 years ago

Starting to code up Euler equation

File size: 4.9 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.