source: trunk/anuga_work/development/Rudy_vandrie_2007/blackbutt/run_blackbutt.py @ 7924

Last change on this file since 7924 was 4441, checked in by ole, 18 years ago

Work with Rudy and Ted 9 May 20007

File size: 5.6 KB
Line 
1"""Script for running a tsunami inundation scenario for Cairns, QLD Australia.
2
3Source data such as elevation and boundary data is assumed to be available in
4directories specified by project.py
5The output sww file is stored in directory named after the scenario, i.e
6slide or fixed_wave.
7
8The scenario is defined by a triangular mesh created from project.polygon,
9the elevation data and a tsunami wave generated by a submarine mass failure.
10
11Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton and
12Nick Bartzis, GA - 2006
13"""
14
15#------------------------------------------------------------------------------
16# Import necessary modules
17#------------------------------------------------------------------------------
18
19# Standard modules
20import os
21import time
22import sys
23
24# Related major packages
25from anuga.shallow_water import Domain
26from anuga.shallow_water import Reflective_boundary
27from anuga.shallow_water import Dirichlet_boundary
28from anuga.shallow_water import Time_boundary
29from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
30from anuga.abstract_2d_finite_volumes.util import file_function
31from anuga.pmesh.mesh_interface import create_mesh_from_regions
32from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf
33from anuga.shallow_water.data_manager import dem2pts
34from anuga.geospatial_data.geospatial_data import Geospatial_data
35from anuga.shallow_water.shallow_water_domain import Inflow
36
37# Application specific imports
38import project                 # Definition of file names and polygons
39
40
41#------------------------------------------------------------------------------
42# Define scenario as either slide or fixed_wave.
43#------------------------------------------------------------------------------
44scenario = 'slide' 
45#scenario = 'fixed_wave'
46
47if os.access(scenario, os.F_OK) == 0:
48    os.mkdir(scenario)
49basename = scenario + 'source'
50
51
52#------------------------------------------------------------------------------
53# Preparation of topographic data
54# Convert ASC 2 DEM 2 PTS using source data and store result in source data
55#------------------------------------------------------------------------------
56
57# Filenames
58basename = 'blackbutt'
59meshname = basename+'.msh'
60
61elevation_file = 'ALL_Grd_Bui_Comb_Smth_Rds_XYZPts_only'
62# Create pts file
63#G = Geospatial_data(elevation_file+'.csv')
64#G.export_points_file(basename+'.pts')
65
66#------------------------------------------------------------------------------
67# Create the triangular mesh based on overall clipping polygon with a tagged
68# boundary and interior regions defined in project.py along with
69# resolutions (maximal area of per triangle) for each polygon
70#------------------------------------------------------------------------------
71
72
73W = 302500
74S = 6171740
75N = 6172270
76E = 304000
77
78bounding_polygon = [[W, S], [E, S], [E, N], [W, N]]
79
80channel_polygon = [[W, S+100], [E, S+100], [E, N-100], [W, N-100]]
81
82interior_regions = [ [channel_polygon, 100] ] # 100 m^2
83create_mesh_from_regions(bounding_polygon,
84                         boundary_tags={'south': [0],
85                                        'east': [1],
86                                        'north': [2],
87                                        'west': [3]},
88                         maximum_triangle_area=400,
89                         interior_regions=interior_regions,
90                         filename=meshname,
91                         use_cache=False,
92                         verbose=True)
93
94
95#------------------------------------------------------------------------------
96# Setup computational domain
97#------------------------------------------------------------------------------
98
99domain = Domain(meshname, use_cache=True, verbose=True)
100
101print domain.statistics()
102domain.set_name(basename) 
103
104
105#------------------------------------------------------------------------------
106# Setup initial conditions
107#------------------------------------------------------------------------------
108
109tide = -100
110domain.set_quantity('stage', tide)
111domain.set_quantity('friction', 0.0, location = 'centroids') 
112domain.set_quantity('elevation', 
113                    filename=basename+'.pts',
114                    use_cache=True,
115                    verbose=True,
116                    alpha=0.01)
117
118
119
120#------------------------------------------------------------------------------
121# Setup specialised forcing terms
122#------------------------------------------------------------------------------
123
124hydrograph = Inflow(center=(300, 300), radius=50,
125             flow=file_function('Q/QPMF_Rot_Sub13.tms'))
126
127domain.forcing_terms.append(hydrograph)
128
129
130#------------------------------------------------------------------------------
131# Setup boundary conditions
132#------------------------------------------------------------------------------
133
134print 'Available boundary tags', domain.get_boundary_tags()
135
136Br = Reflective_boundary(domain)
137Bd = Dirichlet_boundary([tide,0,0])
138
139
140# boundary conditions for slide scenario
141domain.set_boundary({'west': Bd,
142                     'south': Bd,
143                     'east': Transmissive_Momentum_Set_Stage_boundary(domain, lambda t: tide),
144                     'north': Bd})
145   
146
147#------------------------------------------------------------------------------
148# Evolve system through time
149#------------------------------------------------------------------------------
150
151import time
152t0 = time.time()
153
154for t in domain.evolve(yieldstep = 1, finaltime = 3000): 
155    domain.write_time()
156    #domain.write_boundary_statistics(tags='east', quantities='stage')     
157       
158    if t > 2538:
159        hydrograph.flow = 0.0
Note: See TracBrowser for help on using the repository browser.