source: anuga_validation/analytical solutions/Analytical_solution_Yoon_cross_mesh.py @ 5306

Last change on this file since 5306 was 3846, checked in by ole, 17 years ago

Refactored references to domain.filename away.
Use

domain.set_name()
domain.get_name()

instead.

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