source: anuga_work/development/sudi/sw_2d/yoon_parabolic_basin_cross_mesh.py @ 7837

Last change on this file since 7837 was 7837, checked in by mungkasi, 12 years ago

Again, adding some codes for 1d problems on debris avalanche and periodic waves.

File size: 4.5 KB
Line 
1"""Example of shallow water wave equation analytical solution
2consists of a parabolic profile in a parabolic basin. Analytical
3solutiuon to this problem was derived by Thacker
4and used by Yoon and Chou.
5
6   Copyright 2005
7   Christopher Zoppou, Stephen Roberts, ANU, Geoscience Australia
8
9"""
10
11#---------------
12# Module imports
13#from anuga.shallow_water_balanced.swb_domain import Domain, Transmissive_boundary, Reflective_boundary,\
14#     Dirichlet_boundary
15
16from anuga.interface import Domain, Transmissive_boundary, Reflective_boundary,\
17     Dirichlet_boundary
18from math import sqrt, cos, sin, pi
19from anuga.interface import rectangular_cross
20from anuga.utilities.polygon import inside_polygon, is_inside_triangle
21from numpy import asarray
22
23
24#-------------------------------
25# Domain
26n = 20#200
27m = 20#200
28lenx = 8000.0
29leny = 8000.0
30origin = (-4000.0, -4000.0)
31
32points, elements, boundary = rectangular_cross(m, n, lenx, leny, origin)
33domain = Domain(points, elements, boundary)
34
35#----------------
36# Order of scheme
37# Good compromise between
38# limiting and CFL
39#---------------
40domain.set_default_order(2)
41domain.set_timestepping_method(2)
42domain.set_beta(0.7)
43domain.set_CFL(0.6)
44
45domain.smooth = True
46
47#-------------------------------------
48# Provide file name for storing output
49domain.store = False
50domain.format = 'sww'
51domain.set_name('yoon_mesh_second_order_cross')
52print 'Number of triangles = ', len(domain)
53
54#----------------------------------------------------------
55# Decide which quantities are to be stored at each timestep
56domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
57
58#------------------------------------------
59# Reduction operation for get_vertex_values
60from anuga.utilities.numerical_tools import mean
61domain.reduction = mean #domain.reduction = min  #Looks better near steep slopes
62
63#------------------
64# Initial condition
65print 'Initial condition'
66t = 0.0
67D0 = 1.
68L = 2500.
69R0 = 2000.
70g = 9.81
71
72A = (L**4 - R0**4)/(L**4 + R0**4)
73omega = 2./L*sqrt(2.*g*D0)
74T = pi/omega
75
76#------------------
77# Set bed elevation
78def x_slope(x,y):
79    n = x.shape[0]
80    z = 0*x
81    for i in range(n):
82        r = sqrt(x[i]*x[i] + y[i]*y[i])
83        z[i] = -D0*(1.-r*r/L/L)
84    return z
85domain.set_quantity('elevation', x_slope)
86
87#----------------------------
88# Set the initial water level
89def level(x,y):
90    z = x_slope(x,y)
91    n = x.shape[0]
92    h = 0*x
93    for i in range(n):
94        r = sqrt(x[i]*x[i] + y[i]*y[i])
95        h[i] = D0*((sqrt(1-A*A))/(1.-A*cos(omega*t))
96                -1.-r*r/L/L*((1.-A*A)/((1.-A*cos(omega*t))**2)-1.))
97    if h[i] < z[i]:
98        h[i] = z[i]
99    return h
100domain.set_quantity('stage', level)
101
102#---------
103# Boundary
104print 'Boundary conditions'
105R = Reflective_boundary(domain)
106T = Transmissive_boundary(domain)
107D = Dirichlet_boundary([0.0, 0.0, 0.0])
108domain.set_boundary({'left': D, 'right': D, 'top': D, 'bottom': D})
109
110#---------------------------------------------
111# Find triangle that contains the point points
112# and print to file
113
114
115
116points = (0.0, 0.0)
117for n in range(len(domain.triangles)):
118    tri = domain.get_vertex_coordinates(n)
119
120    if is_inside_triangle(points,tri):
121        #print 'Point is within triangle with vertices '+'%s'%tri
122        n_point = n
123
124print 'n_point = ',n_point
125t = domain.triangles[n_point]
126print 't = ', t
127tri = domain.get_vertex_coordinates(n)
128
129filename=domain.get_name()
130file = open(filename,'w')
131
132#----------
133# Evolution
134import time
135t0 = time.time()
136
137time_array = []
138stage_array = []
139Stage     = domain.quantities['stage']
140Xmomentum = domain.quantities['xmomentum']
141Ymomentum = domain.quantities['ymomentum']
142
143for t in domain.evolve(yieldstep = 20.0, finaltime = 17700.0 ):
144    domain.write_time()
145
146    #tri_array = asarray(tri)
147    #t_array = asarray([[0,1,2]])
148    #interp = Interpolation(tri_array,t_array,[points])
149
150
151    stage     = Stage.get_values(location='centroids',indices=[n_point])
152    xmomentum = Xmomentum.get_values(location='centroids',indices=[n_point])
153    ymomentum = Ymomentum.get_values(location='centroids',indices=[n_point])
154    #print '%10.6f   %10.6f  %10.6f   %10.6f\n'%(t,stage,xmomentum,ymomentum)
155    file.write( '%10.6f   %10.6f  %10.6f   %10.6f\n'%(t,stage,xmomentum,ymomentum) )
156
157    time_array.append(t)
158    stage_array.append(stage)
159
160file.close()
161print 'That took %.2f seconds' %(time.time()-t0)
162
163
164from pylab import *
165ion()
166hold(False)
167plot(time_array, stage_array, 'r.-')
168#title('Gauge %s' %name)
169xlabel('time(s)')
170ylabel('stage (m)')
171#legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
172#savefig(name, dpi = 300)
173
174raw_input('Next')
175show()
176
177
Note: See TracBrowser for help on using the repository browser.