source: trunk/anuga_work/development/sudi/sw_2d/yoon_parabolic_basin_cross_mesh.py @ 7886

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

Copied sudi's directory

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