source: trunk/anuga_validation/validation_tests/parabolic_basin.py @ 8102

Last change on this file since 8102 was 8102, checked in by steve, 13 years ago

Adding in files

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