source: anuga_validation/analytical_solutions/macdonald_import_mesh.py @ 7636

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

Cleaning up directory

File size: 3.1 KB
Line 
1"""Example of shallow water wave equation analytical solution
2consists of a symmetrical converging channel with friction and bed slope.
3
4Specific methods pertaining to the 2D shallow water equation
5are imported from shallow_water
6for use with the generic finite volume framework
7
8   Copyright 2004
9   Christopher Zoppou, Stephen Roberts, Ole Nielsen, Duncan Gray
10   Geoscience Australia, ANU
11   
12Specific methods pertaining to the 2D shallow water equation
13are imported from shallow_water
14for use with the generic finite volume framework
15
16Conserved quantities are h, uh and vh stored as elements 0, 1 and 2 in the
17numerical vector named conserved_quantities.
18"""
19
20#-----------------------------
21# Module imports
22import sys, string
23from os import sep
24sys.path.append('..'+sep+'pyvolution')
25
26from shallow_water import Transmissive_boundary, Reflective_boundary, \
27     Dirichlet_boundary
28from shallow_water import Domain
29from pmesh2domain import pmesh_to_domain_instance
30from Numeric import zeros, Float
31
32
33#------------------------------
34# Domain
35filename = 'MacDonald_77541.tsh'
36print 'Creating domain from', filename
37domain = pmesh_to_domain_instance(filename, Domain)
38print 'Number of triangles = ', len(domain)
39
40domain.default_order = 1
41domain.smooth = True
42
43#-------------------------------------
44# Provide file name for storing output
45domain.store = True     #Store for visualisation purposes
46domain.format = 'sww'   #Native netcdf visualisation format
47domain.set_name('MacDonald_77541')
48
49#-------------
50#Bed Elevation
51fid = open('MacDonald_bed.xya')
52lines = fid.readlines()
53n_bed = len(lines)
54z_bed = zeros(n_bed, Float)
55x_bed = zeros(n_bed, Float)
56
57for line in range(n_bed):
58    value = string.split(lines[line])
59    x_bed[line] = float(value[0])
60    z_bed[line] = float(value[1])
61   
62#-----------------
63#Set bed-elevation
64def x_slope(x,y):
65    n = x.shape[0]
66    z = 0*x
67    for i in range(n):
68        for j in range(n_bed-1):
69            if x[i] >= x_bed[j] and x[i] < x_bed[j+1]:
70                z[i] = z_bed[j] + (z_bed[j+1] - z_bed[j])/(x_bed[j+1] - x_bed[j])*(x[i] - x_bed[j])
71            return z
72
73domain.set_quantity('elevation', x_slope)
74   
75#---------------------------------------------------------
76#Decide which quantities are to be stored at each timestep
77domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
78
79#Reduction operation for get_vertex_values             
80from anuga.pyvolution.util import mean
81domain.reduction = mean
82
83#--------------------
84# Boundary Conditions
85upstream_depth = 5.035 - 4.393
86downstream_depth = 1.5 - 0
87discharge = 20
88tags = {}
89tags['Upstream'] = Dirichlet_boundary([upstream_depth, 2, 0.0])
90tags['external'] = Reflective_boundary(domain) 
91tags['Downstream'] = Dirichlet_boundary([downstream_depth, 2, 0.0])
92domain.set_boundary(tags)
93
94#---------
95# friction
96domain.set_quantity('friction', 0.02)
97
98#-----------------
99#Initial condition
100domain.set_quantity('elevation', 0.0)
101domain.set_quantity('stage', 0.2)
102   
103#-----------------
104#Evolution
105import time
106t0 = time.time()
107for t in domain.evolve(yieldstep = 5, finaltime = 500):
108    domain.write_time()
109   
110print 'That took %.2f seconds' %(time.time()-t0)
111
112   
Note: See TracBrowser for help on using the repository browser.