source: anuga_work/development/analytical_solutions/Analytical_solution_Sampson_import_mesh.py @ 5029

Last change on this file since 5029 was 5029, checked in by steve, 16 years ago

Fixing up directory name

File size: 2.9 KB
Line 
1"""Example of shallow water wave equation analytical solution
2consists of a flat water surface profile in a parabolic basin
3with linear friction. The analytical solution was derived by Sampson in 2002.
4
5   Copyright 2004
6   Christopher Zoppou, Stephen Roberts, Ole Nielsen, Duncan Gray
7   Geoscience Australia
8   
9Specific methods pertaining to the 2D shallow water equation
10are imported from shallow_water
11for use with the generic finite volume framework
12
13Conserved quantities are h, uh and vh stored as elements 0, 1 and 2 in the
14numerical vector named conserved_quantities.
15"""
16
17######################
18# Module imports
19#
20
21import sys
22from os import sep
23sys.path.append('..'+sep+'pyvolution')
24
25from shallow_water import Domain, Dirichlet_boundary, gravity, linear_friction
26from math import sqrt, cos, sin, pi, exp
27from mesh_factory import strang_mesh
28from quantity import Quantity
29
30
31######################
32# Domain
33# Strang_domain will search through the file and test to see if there are
34# two or three entries. Two entries are for points and three for triangles.
35
36
37points, elements = strang_mesh('yoon_circle.pt')
38domain = Domain(points, elements) 
39
40domain.default_order = 2
41domain.smooth = True
42
43domain.quantities['linear_friction'] = Quantity(domain)
44
45#Reconstruct forcing terms with linear friction instead og manning friction
46domain.forcing_terms = []
47domain.forcing_terms.append(gravity)
48domain.forcing_terms.append(linear_friction)
49
50print domain.forcing_terms
51
52# Provide file name for storing output
53domain.store = True
54domain.format = 'sww'
55domain.filename = 'sampson_strang_second_order'
56
57print 'Number of triangles = ', len(domain)
58
59
60#Reduction operation for get_vertex_values             
61from anuga.pyvolution.util import mean
62domain.reduction = mean
63#domain.reduction = min  #Looks better near steep slopes
64
65
66######################
67#Initial condition
68#
69print 'Initial condition'
70t = 0.0
71h0 = 10.
72a = 3000.
73g = 9.81
74tau =0.001
75B = 5
76A = 0
77p = sqrt(8*g*h0/a/a)
78s = sqrt(p*p-tau*tau)/2
79t = 0.
80
81
82#Set bed-elevation and friction
83def x_slope(x,y):
84    n = x.shape[0]
85    z = 0*x
86    for i in range(n):
87        z[i] = h0 - h0*(1.0 -x[i]*x[i]/a/a - y[i]*y[i]/a/a)
88    return z
89
90domain.set_quantity('elevation', x_slope)
91domain.set_quantity('linear_friction', tau)
92
93
94#Set the water stage
95def stage(x,y):
96    z = x_slope(x,y)
97    n = x.shape[0]
98    h = 0*x
99    for i in range(n):
100        h[i] = h0-B*B*exp(-tau*t)/2/g-1/g*(exp(-tau*t/2)*(B*s*cos(s*t) \
101               +tau*B/2*sin(s*t)))*x[i] \
102          -1/g*(exp(-tau*t/2)*(B*s*sin(s*t)-tau*B/2*cos(s*t)))*y[i]         
103    if h[i] < z[i]:
104        h[i] = z[i]
105    return h   
106
107domain.set_quantity('stage', stage)
108
109
110############
111#Boundary
112domain.set_boundary({'exterior': Dirichlet_boundary([0.0, 0.0, 0.0])})
113   
114
115######################
116#Evolution
117import time
118t0 = time.time()
119for t in domain.evolve(yieldstep = 10.0, finaltime = 5000):
120    domain.write_time()
121
122print 'That took %.2f seconds' %(time.time()-t0)
123
124   
Note: See TracBrowser for help on using the repository browser.