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

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

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

File size: 5.0 KB
Line 
1"""Example of shallow water wave equation analytical solution of the
2circular hydraulic jump experimental data treated as a two-dimensional solution.
3
4   Copyright 2005
5   Christopher Zoppou, Stephen Roberts
6   Geoscience Australia, ANU
7"""
8
9#-------------------------------
10# Set up path and module imports
11import sys
12from os import sep
13sys.path.append('..'+sep+'pyvolution')
14
15from shallow_water import Domain, Dirichlet_Discharge_boundary
16from shallow_water import Transmissive_Momentum_Set_Stage_boundary, Dirichlet_boundary
17from math import pi, sqrt
18from mesh_factory import strang_mesh
19
20
21#---------
22# Geometry
23bed = ([.519, .519, .519, .519, .5192, .5194, .5196, .520, .5207, .5215, .5233, .5233])
24distance = ([.08, .10, .11, .16, .21, .26, .31, .36, .41, .46, .50, .52])
25n_bed = 12
26
27#---------
28# Case A.4
29Q = 9.985/1000.0
30wh0 = Q/(2.0*pi*0.1)
31stage0 = bed[2] + 0.005
32wh1 = -Q/(2.0*pi*0.503)
33stage1 = 0.562
34Manning = 0.009
35
36#------------------
37# Set up the domain
38# Strang_domain will search through the file and test to see if there are
39# two or three entries. Two entries are for points and three for triangles.
40points, elements = strang_mesh('circular.pt')
41domain = Domain(points, elements)
42
43print "Number of triangles = ", len(domain)
44
45#----------------------
46# Set a default tagging
47
48
49for id, face in domain.boundary:
50    domain.boundary[(id,face)] = 'outer'
51    point = domain.get_vertex_coordinate(id,(face+1)%3)
52    radius2 = point[0]*point[0] + point[1]*point[1]
53    typical_outer = (id,face)
54    if radius2 < 0.1:
55        domain.boundary[(id,face)] = 'inner'
56        typical_inner = (id,face)
57
58
59#-------------------------------------
60# Provide file name for storing output
61#domain.visualise = True
62domain.store = True
63domain.format = 'sww'
64domain.filename = 'circular_second_order'
65
66#------------------------------------------
67# Reduction operation for get_vertex_values
68from anuga.pyvolution.util import mean
69domain.reduction = mean
70#domain.reduction = min  #Looks better near steep slopes
71
72#---------------------------
73# Function for bed-elevation
74def bed_z(x,y):
75    n = x.shape[0]
76    z = 0*x
77    for i in range(n):
78        r = sqrt(x[i]*x[i]+y[i]*y[i])
79        for j in range(n_bed-1):
80            if distance[j] <= r:
81                if distance[j+1] > r:
82                    z[i] = bed[j] + (bed[j+1] - bed[j])/(distance[j+1] - distance[j])*(r - distance[j])
83    return z
84
85domain.set_quantity('elevation', bed_z)
86
87#---------
88# Friction
89domain.set_quantity('friction', Manning)
90
91#---------------------------------
92# Function for initial water elevation
93# (stage)
94def level(x,y):
95    z = bed_z(x,y)
96    n = x.shape[0]
97    stage = 0*x
98    for i in range(n):
99        stage[i] = stage0
100    return stage
101
102
103#def outflow_stage(t):
104#    return [stage1, 0 , 0]
105
106
107domain.set_quantity('stage', level)
108
109#---------------------------
110# Set up boundary conditions
111DD_BC_INNER = Dirichlet_Discharge_boundary(domain, stage0, wh0)
112DD_BC_OUTER = Dirichlet_Discharge_boundary(domain, stage1, wh1)
113
114domain.set_boundary({'inner': DD_BC_INNER, 'outer': DD_BC_OUTER})
115
116#------------------
117# Order of accuracy
118domain.default_order = 1
119domain.CFL = 0.75
120#domain.beta_w = 0.5
121#domain.beta_h = 0.2
122domain.smooth = True
123
124
125# domain.initialise_visualiser()
126# domain.visualiser.coloring['stage'] = True
127# domain.visualiser.scale_z['stage'] = 2.0
128# domain.visualiser.scale_z['elevation'] = 0.05
129
130
131#domain.initialise_visualiser()
132#    #domain.visualiser.coloring['stage'] = True
133#domain.visualiser.scale_z['stage'] = 2.0
134#domain.visualiser.scale_z['elevation'] = 0.05
135#
136#
137#from realtime_visualisation_new import Visualiser
138##vxmom = Visualiser(domain,title='xmomentum',scale_z=10.0)
139##vymom = Visualiser(domain,title='ymomentum',scale_z=10.0)
140
141stage = domain.quantities['stage']
142
143#----------
144# Evolution
145import time
146
147f = open("circular_hydraulic_jump_true.txt","w")
148
149t0 = time.time()
150for t in domain.evolve(yieldstep = .1, finaltime = 3.0):
151    domain.write_time()
152
153    exp = '(xmomentum**2 + ymomentum**2)**0.5'
154    radial_momentum = domain.create_quantity_from_expression(exp)
155
156    print 'outer stage      ', stage.get_values(location='vertices',
157                                            indices=[typical_outer[0]])
158    print '      radial mom ', \
159          radial_momentum.get_values(location='centroids',
160                                     indices=[typical_outer[0]])
161
162    print 'inner stage      ', stage.get_values(location='centroids',
163                                            indices=[typical_inner[0]])
164    print '      radial mom ', \
165          radial_momentum.get_values(location='centroids',
166                                     indices=[typical_inner[0]])
167
168#    f.write('time = %25.15e wall clock time %g \n' % (domain.time, time.time()))
169    f.write('%10.3f %25.15e %25.15e %25.15e %25.15e \n' % (domain.time, inner_stage, inner_radial_mom, outer_stage, outer_radial_mom))
170
171    f.write('time = %25.15e wall clock time %g \n' % (domain.time, time.time()))
172    f.write('%g \n' % stage.get_values(location='centroids',
173                                   indices=[typical_outer[0]])[0])
174
175f.close()
176
177print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.