source: development/analytical solutions/Analytical_solution_Yoon_import_mesh.py @ 3514

Last change on this file since 3514 was 3514, checked in by duncan, 18 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.0 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, Ole Nielsen
8   ANU, Geoscience Australia
9   
10"""
11
12#---------------
13# Module imports
14import sys
15from os import sep
16sys.path.append('..'+sep+'pyvolution')
17from shallow_water import Domain, Dirichlet_boundary
18from math import sqrt, cos, sin, pi
19from mesh_factory import strang_mesh
20from anuga.pyvolution.util import inside_polygon
21from Numeric import asarray
22from least_squares import Interpolation
23
24#-------------------------------
25# Set up the domain of triangles
26# Strang_domain will search
27# through the file and test to
28# see if there are two or three
29# entries. Two entries are for
30# points and three for triangles.
31points, elements = strang_mesh('yoon_circle.pt')
32domain = Domain(points, elements) 
33
34#----------------
35# Order of scheme
36domain.default_order = 2
37
38domain.smooth = True
39
40#-------------------------------------
41# Provide file name for storing output
42domain.store = False
43domain.format = 'sww'
44domain.filename = 'yoon_mesh_second_order.2'
45print 'Number of triangles = ', len(domain)
46
47#----------------------------------------------------------
48# Decide which quantities are to be stored at each timestep
49domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
50
51#------------------------------------------
52# Reduction operation for get_vertex_values
53from anuga.pyvolution.util import mean
54domain.reduction = mean #domain.reduction = min  #Looks better near steep slopes
55
56
57#------------------
58# Initial condition
59print 'Initial condition'
60t = 0.0
61D0 = 1.
62L = 2500.
63R0 = 2000.
64g = 9.81
65
66A = (L**4 - R0**4)/(L**4 + R0**4)
67omega = 2./L*sqrt(2.*g*D0)
68T = pi/omega
69
70#------------------
71# Set bed elevation
72def x_slope(x,y):
73    n = x.shape[0]
74    z = 0*x
75    for i in range(n):
76        r = sqrt(x[i]*x[i] + y[i]*y[i])
77        z[i] = -D0*(1.-r*r/L/L)
78    return z
79domain.set_quantity('elevation', x_slope)
80
81#----------------------------
82# Set the initial water level
83def level(x,y):
84    z = x_slope(x,y)
85    n = x.shape[0]
86    h = 0*x
87    for i in range(n):
88        r = sqrt(x[i]*x[i] + y[i]*y[i])
89        h[i] = D0*((sqrt(1-A*A))/(1.-A*cos(omega*t))
90                -1.-r*r/L/L*((1.-A*A)/((1.-A*cos(omega*t))**2)-1.))
91    if h[i] < z[i]:
92        h[i] = z[i]
93    return h   
94domain.set_quantity('stage', level)
95
96
97#---------
98# Boundary
99print 'Boundary conditions'
100domain.set_boundary({'exterior': Dirichlet_boundary([0.0, 0.0, 0.0])})
101
102#---------------------------------------------
103# Find triangle that contains the point points
104# and print to file
105points = [0.,0.]
106for n in range(len(domain.triangles)):
107    t = domain.triangles[n]
108    tri = [domain.coordinates[t[0]],domain.coordinates[t[1]],domain.coordinates[t[2]]]
109
110    if inside_polygon(points,tri):
111        print 'Point is within triangle with vertices '+'%s'%tri
112        n_point = n
113t = domain.triangles[n_point]
114tri = [domain.coordinates[t[0]],domain.coordinates[t[1]],domain.coordinates[t[2]]]
115
116filename=domain.filename
117file = open(filename,'w')
118   
119#----------
120# Evolution
121import time
122t0 = time.time()
123for t in domain.evolve(yieldstep = 20.0, finaltime = 3000 ):
124    domain.write_time()
125
126    tri_array = asarray(tri)
127    t_array = asarray([[0,1,2]])
128    interp = Interpolation(tri_array,t_array,[points])
129
130    stage = domain.get_quantity('stage').get_values()[n_point]
131    xmomentum = domain.get_quantity('xmomentum').get_values()[n_point]
132    ymomentum = domain.get_quantity('ymomentum').get_values()[n_point]
133
134    interp_stage = interp.interpolate([[stage[0]],[stage[1]],[stage[2]]])
135    interp_xmomentum = interp.interpolate([[xmomentum[0]],[xmomentum[1]],[xmomentum[2]]])
136    interp_ymomentum = interp.interpolate([[ymomentum[0]],[ymomentum[1]],[ymomentum[2]]])
137   
138    file.write( '%10.6f   %10.6f  %10.6f   %10.6f\n'%(t,interp_stage[0][0],interp_xmomentum[0][0],interp_ymomentum[0][0]) )
139
140file.close()   
141print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.