source: inundation/analytical solutions/Analytical_solution_circular_hydraulic_jump.py @ 2001

Last change on this file since 2001 was 2001, checked in by chris, 19 years ago
File size: 4.6 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)
31h0 = bed[2] + 0.005
32wh1 = -Q/(2.0*pi*0.5)
33h1 = 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 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 depth
93def level(x,y):
94    z = bed_z(x,y)
95    n = x.shape[0]
96    w = 0*x
97    for i in range(n):
98        w[i] = h0
99        h = w[i] - z[i]
100    return w
101
102
103def outflow_height(t):
104    return [h1, 0 , 0]
105
106
107domain.set_quantity('stage', level)
108
109#---------------------------
110# Set up boundary conditions
111DD_BC_INNER = Dirichlet_Discharge_boundary(domain, h0, wh0)
112DD_BC_OUTER = Dirichlet_Discharge_boundary(domain, h1, 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# from realtime_visualisation_new import Visualiser
132# vxmom = Visualiser(domain,title='xmomentum',scale_z=10.0)
133# vymom = Visualiser(domain,title='ymomentum',scale_z=10.0)
134
135w = domain.quantities['stage']
136
137#----------
138# Evolution
139import time
140
141f = open("circular_hydraulic_jump_true.txt","w")
142
143t0 = time.time()
144for t in domain.evolve(yieldstep = 1.0, finaltime = 1000.):
145    domain.write_time()
146
147    exp = '(xmomentum**2 + ymomentum**2)**0.5'
148    radial_momentum = domain.create_quantity_from_expression(exp)
149
150    outer_stage = w.get_values(location='centroids', indices=[typical_outer[0]])
151    outer_radial_mom = radial_momentum.get_values(location='centroids', indices=[typical_outer[0]])
152   
153    inner_stage = w.get_values(location='centroids', indices=[typical_inner[0]])
154    inner_radial_mom = radial_momentum.get_values(location='centroids', indices=[typical_inner[0]])
155
156    print inner_stage
157    print inner_radial_mom
158    print outer_stage
159    print outer_radial_mom
160
161#    f.write('time = %25.15e wall clock time %g \n' % (domain.time, time.time()))
162    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))
163
164
165f.close()
166
167print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.