source: anuga_work/development/analytical solutions/Analytical_solution_Yoon_cross_mesh.py @ 3970

Last change on this file since 3970 was 3514, checked in by duncan, 19 years ago

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

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.filename = '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.filename
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.