source: anuga_validation/analytical solutions/Analytical solution_oblique_shock_01.py @ 5306

Last change on this file since 5306 was 3846, checked in by ole, 17 years ago

Refactored references to domain.filename away.
Use

domain.set_name()
domain.get_name()

instead.

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
26from anuga.shallow_water import Domain
27from anuga.shallow_water import Transmissive_boundary, Reflective_boundary,\
28     Dirichlet_boundary , Transmissive_Momentum_Set_Stage_boundary
29
30from math import sqrt, cos, sin, pi
31from mesh_factory import oblique_cross
32
33
34######################
35# Domain
36#
37leny = 30.
38lenx = 40.
39f = 2
40n = 50*f
41m = 80*f
42theta = 25
43h_bc = 1.5
44p_bc = 5.0
45
46points, elements, boundary = oblique_cross(m, n, lenx, leny, theta = theta)
47domain = Domain(points, elements, boundary)
48
49
50print 'Number of Elements = ',domain.number_of_elements
51# Order of solver
52domain.default_order=2
53domain.beta_w = 0.8
54
55# Store output
56domain.store=True
57
58# Provide file name for storing output
59domain.set_name('oblique_cross_%g_%g_%g_%g_%g' % (n,m,theta,h_bc,p_bc))
60print domain.get_name()
61
62# Output format
63domain.format="sww" #NET.CDF binary format
64                    # "dat" for ASCII
65
66
67# Visualization smoothing
68domain.smooth=True
69domain.visualise=False
70
71#######################
72#Bed-slope and friction
73class ConstantFunctionT:
74
75    def __init__(self,value):
76        self.value = value
77
78    def __call__(self,t):
79        return self.value
80
81
82shock_hb = ConstantFunctionT(h_bc)
83shock_hh = ConstantFunctionT(h_bc)
84shock_pb = ConstantFunctionT(p_bc)
85shock_pp = ConstantFunctionT(p_bc)
86
87
88def x_slope(x, y):
89    return 0*x
90
91
92
93def shock_h(x,y):
94    n = x.shape[0]
95    w = 0*x
96    for i in range(n):
97        if x[i]<0:
98            w[i] = shock_hh(0.0)
99        else:
100            w[i] = 1.0
101    return w
102
103
104def shock_p(x,y):
105    n = x.shape[0]
106    w = 0*x
107    for i in range(n):
108        if x[i]<0:
109            w[i] = shock_pp(0.0)
110        else:
111            w[i] = 0.0
112    return w
113
114
115
116
117
118domain.set_quantity('elevation', x_slope)
119domain.set_quantity('friction', 0.0)
120
121######################
122# Boundary conditions
123#
124R = Reflective_boundary(domain)
125T = Transmissive_boundary(domain)
126D = Dirichlet_boundary([shock_hb(0.0) , shock_pb(0.0), 0.0])
127Bts = Transmissive_Momentum_Set_Stage_boundary(domain, shock_hb)
128
129domain.set_boundary({'left': D, 'right': T, 'top': R, 'bottom': R})
130
131######################
132#Initial condition
133
134domain.set_quantity('stage', shock_h )
135domain.set_quantity('xmomentum',shock_p )
136
137class SetValueWhere:
138
139    def __init__(self,quantity,value=0,xrange=0.0):
140        self.quantity =  quantity
141        self.value = value
142        self.xrange = xrange
143
144    def __call__(self,x,y):
145        n = x.shape[0]
146        w = self.quantity.centroid_values
147        for i in range(n):
148            if x[i]<self.xrange:
149                w[i] = self.value
150        return w
151
152
153Stage = domain.quantities['stage']
154Xmom =  domain.quantities['xmomentum']
155Ymom =  domain.quantities['ymomentum']
156
157#for id, face in domain.boundary:
158#    print '(',id,',',face,') ',domain.boundary[(id,face)]
159
160######################
161#Evolution
162import time
163t0 = time.time()
164for t in domain.evolve(yieldstep = 0.1, finaltime = 20.0):
165    domain.write_time()
166    id = 3399
167    print Stage.get_values(location='centroids',indices=[id])
168    print Xmom.get_values(location='centroids',indices=[id])
169    print Ymom.get_values(location='centroids',indices=[id])
170    vstage =  Stage.get_values(location='centroids',indices=[id])
171    vxmom  =  Xmom.get_values(location='centroids',indices=[id])
172    id = 12719
173    print Stage.get_values(location='centroids',indices=[id])
174    print Xmom.get_values(location='centroids',indices=[id])
175    print Ymom.get_values(location='centroids',indices=[id])
176    tclean = 2.0
177#    if t > tclean-0.09 and t < tclean+0.01:
178#        print 'Cleaning up Initial profile'
179#        setstage = SetValueWhere(quantity=Stage,value=vstage[0],xrange=10)
180#        setxmom  = SetValueWhere(quantity=Xmom,value=vxmom[0],xrange=10)
181#        Stage.set_values(setstage,location='centroids')
182#        Xmom.set_values(setxmom,location='centroids')
183#        Dnew = Dirichlet_boundary([vstage[0] , vxmom[0], 0.0])
184#        domain.set_boundary({'left': Dnew, 'right': T, 'top': R, 'bottom': R})
185
186print 'That took %.2f seconds' %(time.time()-t0)
187
188#FIXME: Compute average water depth on either side of shock and compare
189#to expected values. And also Froude numbers.
190
191
192#print "saving file?"
193#im = ImageGrab.grab()
194#im.save("ccube.eps")
Note: See TracBrowser for help on using the repository browser.