source: production/wollongong_2006/run_flagstaff.py @ 3046

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

changed parameters

File size: 4.6 KB
Line 
1"""Script for running a tsunami inundation scenario for Flagstaff pt,
2Wollongong harbour, NSW, Australia.
3
4Source data such as elevation and boundary data is assumed to be available in
5directories specified by project.py
6
7The scenario is defined by a triangular mesh created from project.polygon,
8the elevation data and a hypothetical boundary condition.
9
10Ole Nielsen and Duncan Gray, GA - 2005, Nick Bartzis and Jane Sexton, GA - 2006
11"""
12
13
14#------------------------------------------------------------------------------
15# Import necessary modules
16#------------------------------------------------------------------------------
17
18
19# Standard modules
20import os, sys, time 
21from os import sep
22from os.path import dirname, basename
23
24# Related major packages
25from pyvolution.shallow_water import Domain
26from pyvolution.shallow_water import Dirichlet_boundary
27from pyvolution.shallow_water import Time_boundary
28from pyvolution.shallow_water import Reflective_boundary
29from pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts
30from pmesh.mesh_interface import create_mesh_from_regions
31from pmesh.mesh import importUngenerateFile, Segment
32
33# Application specific imports
34import project
35
36
37#------------------------------------------------------------------------------
38# Preparation of topographic data
39#
40# Convert ASC 2 DEM 2 PTS using source data and store result in source data
41#------------------------------------------------------------------------------
42
43# Create DEM from asc data
44convert_dem_from_ascii2netcdf(project.demname, use_cache=True, verbose=True)
45
46# Create pts file from DEM
47dem2pts(project.demname,
48        easting_min=project.xllcorner,
49        easting_max=project.xurcorner,
50        northing_min=project.yllcorner,
51        northing_max= project.yurcorner,
52        use_cache=True, 
53        verbose=True)
54
55
56#------------------------------------------------------------------------------
57# Create the triangular mesh based on overall clipping polygon with a tagged
58# boundary and interior regions defined in project.py along with
59# resolutions (maximal area of per triangle) for each polygon
60#------------------------------------------------------------------------------
61
62# Generate basic mesh
63mesh = create_mesh_from_regions(project.bounding_polygon,
64                                boundary_tags=project.boundary_tags,
65                                maximum_triangle_area=project.base_resolution,
66                                interior_regions=project.interior_regions)
67
68# Add buildings
69dict = importUngenerateFile(project.buildings_filename)
70Segment.set_default_tag('wall') # This should bind to a Reflective boundary
71mesh.addVertsSegs(dict)
72
73# Generate and write mesh to file
74mesh.generate_mesh(maximum_triangle_area=project.base_resolution,
75                   verbose=True)
76
77mesh.export_mesh_file(project.mesh_filename)
78
79
80#------------------------------------------------------------------------------
81# Setup computational domain
82#------------------------------------------------------------------------------
83
84domain = Domain(project.mesh_filename, use_cache = False, verbose = True)
85print domain.statistics()
86
87domain.set_name(project.basename)
88domain.set_datadir(project.outputdir)
89domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
90
91
92#------------------------------------------------------------------------------
93# Setup initial conditions
94#------------------------------------------------------------------------------
95
96domain.set_quantity('stage', project.initial_sealevel)
97domain.set_quantity('friction', 0.03)
98domain.set_quantity('elevation', 
99                    filename=project.demname + '.pts',
100                    use_cache=True,
101                    verbose=True)
102
103
104#------------------------------------------------------------------------------
105# Setup boundary conditions
106#------------------------------------------------------------------------------
107
108D = Dirichlet_boundary([project.initial_sealevel, 0, 0])
109R = Reflective_boundary(domain)
110W = Time_boundary(domain = domain,
111                  f=lambda t: [project.initial_sealevel + (60<t<480)*6, 0, 0])
112
113domain.set_boundary({'exterior': D,
114                     'side': D,
115                     'wall': R,
116                     'ocean': W})
117
118
119                             
120#------------------------------------------------------------------------------
121# Evolve system through time
122#------------------------------------------------------------------------------
123
124t0 = time.time()
125for t in domain.evolve(yieldstep = 1, finaltime = 1200): 
126    domain.write_time()
127    domain.write_boundary_statistics(tags = 'ocean')     
128 
129print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.