source: anuga_validation/analytical solutions/Analytical_solution_MacDonald_import_mesh_wall0.py @ 3846

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

Refactored references to domain.filename away.
Use

domain.set_name()
domain.get_name()

instead.

File size: 3.6 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. Existance of 'MacDonald_77541.tsh' assumed.
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, Add_value_to_region
29from pmesh2domain import pmesh_to_domain_instance
30from Numeric import zeros, Float
31
32#-------
33# Domain
34filename = 'MacDonald_91328_wall.tsh'
35print 'Creating domain from', filename
36domain = pmesh_to_domain_instance(filename, Domain)
37print 'Number of triangles = ', len(domain)
38
39#-------------------------------------
40# Provide file name for storing output
41domain.store = True     #Store for visualisation purposes
42domain.format = 'sww'   #Native netcdf visualisation format
43domain.set_name('MacDonald_first_order_wall')
44
45#--------------
46# Bed Elevation
47fid = open('MacDonald_bed.xya')
48lines = fid.readlines()
49n_bed = len(lines)
50z_bed = zeros(n_bed, Float)
51x_bed = zeros(n_bed, Float)
52
53for line in range(n_bed):
54    value = string.split(lines[line])
55    x_bed[line] = float(value[0])
56    z_bed[line] = float(value[1])
57   
58#------------------
59# Set bed elevation
60def bed_elevation(x,y):
61    n = x.shape[0]
62    z = 0*x
63   
64    for i in range(n):
65#       print i
66        for j in range(n_bed-1):
67            if x[i] >= x_bed[j] and x[i] < x_bed[j+1]:
68                z[i] = z_bed[j] + (z_bed[j+1] - z_bed[j])/(x_bed[j+1] - x_bed[j])*(x[i] - x_bed[j])
69    return z
70
71# Bed elevation         
72domain.set_quantity('elevation', bed_elevation)
73
74#-----------------------------------------
75# Set initial water surface stage which is
76# bed elevation plus an arbitrary 0.2   
77#---------------------------
78# Add value to tagged region
79domain.set_region(Add_value_to_region('bed', 'stage', 0.2,location='unique vertices', initial_quantity='elevation'))
80
81#---------------------------
82# Add value to tagged region
83domain.set_region(Add_value_to_region('floodplain', 'elevation', 2,location='unique vertices', initial_quantity='elevation'))
84   
85#----------------------------------------------------------
86# Decide which quantities are to be stored at each timestep
87domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
88
89#------------------------------------------
90# Reduction operation for get_vertex_values             
91from anuga.pyvolution.util import mean
92domain.reduction = mean
93
94#--------------------
95# Boundary Conditions
96upstream_depth = 5.035 
97downstream_depth = 1.5
98discharge = 20
99tags = {}
100tags['upstream'] = Dirichlet_boundary([upstream_depth, 2, 0.0])
101tags['exterior'] = Reflective_boundary(domain) 
102tags['downstream'] = Dirichlet_boundary([downstream_depth, 2, 0.0])
103domain.set_boundary(tags)
104
105#---------
106# Friction
107domain.set_quantity('friction', 0.02)
108
109# -------------------
110# Set order of scheme
111domain.default_order = 1
112domain.smooth = True
113
114#----------
115# Evolution
116import time
117t0 = time.time()
118for t in domain.evolve(yieldstep = 100, finaltime = 1500.):
119    domain.write_time()
120   
121print 'That took %.2f seconds' %(time.time()-t0)
122
123   
Note: See TracBrowser for help on using the repository browser.