source: anuga_validation/analytical solutions/circular_hydraulic_jump.py @ 7626

Last change on this file since 7626 was 7626, checked in by steve, 14 years ago

Cleaning up directory

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