source: branches/numpy_anuga_validation/analytical solutions/Analytical_solution_MacDonald_import_mesh_wall.py @ 7080

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

Refactored references to domain.filename away.
Use

domain.set_name()
domain.get_name()

instead.

File size: 3.3 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, Add_value_to_region
29from pmesh2domain import pmesh_to_domain_instance
30from Numeric import zeros, Float
31
32
33#------------------------------
34# Domain
35filename = 'MacDonald_77541_wall.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_first_order_wall')
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# Add value to taged region
77domain.set_region(Add_value_to_region('floodplain', 'elevation', 2, location='unique vertices', initial_quantity='elevation'))
78
79#---------------------------------------------------------
80#Decide which quantities are to be stored at each timestep
81domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
82
83#Reduction operation for get_vertex_values             
84from anuga.pyvolution.util import mean
85domain.reduction = mean
86
87#--------------------
88# Boundary Conditions
89upstream_depth = 5.035 - 4.393
90downstream_depth = 1.5
91discharge = 20
92tags = {}
93tags['Upstream'] = Dirichlet_boundary([upstream_depth, 2.0, 0.0])
94tags['Downstream'] = Dirichlet_boundary([downstream_depth, 0.0, 0.0])
95tags['exterior'] = Reflective_boundary(domain) 
96domain.set_boundary(tags)
97
98#---------
99# friction
100domain.set_quantity('friction', 0.02)
101
102#-----------------
103#Initial condition
104#domain.set_quantity('elevation', 0.0)
105domain.set_quantity('stage', 0.2)
106   
107#-----------------
108#Evolution
109import time
110t0 = time.time()
111for t in domain.evolve(yieldstep = .1, finaltime = 5):
112    domain.write_time()
113   
114print 'That took %.2f seconds' %(time.time()-t0)
115
116   
Note: See TracBrowser for help on using the repository browser.