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

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

Cleaning up directory

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