source: anuga_work/production/busselton/standardised_version/run_busselton.py @ 6259

Last change on this file since 6259 was 6259, checked in by ole, 10 years ago

Work on simplifying model script

File size: 6.0 KB
Line 
1"""Script for running a tsunami inundation scenario for busselton, WA, Australia.
2
3The scenario is defined by a triangular mesh created from project.polygon,
4the elevation data is compiled into a pts file through build_busselton.py
5and a simulated tsunami is generated through an sts file from build_boundary.py.
6
7Input: sts file (build_boundary.py for respective event)
8       pts file (build_busselton.py)
9       information from project file
10Outputs: sww file stored in project.output_run_time_dir
11The export_results_all.py and get_timeseries.py is reliant
12on the outputs of this script
13
14Ole Nielsen and Duncan Gray, GA - 2005, Jane Sexton, Nick Bartzis, GA - 2006
15Ole Nielsen, Jane Sexton and Kristy Van Putten - 2008
16"""
17
18# Note repeated use of domain id
19# Note naming of fundamental datasets and directories
20
21
22#------------------------------------------------------------------------------
23# Import necessary modules
24#------------------------------------------------------------------------------
25
26# Standard modules
27from os import sep
28import os
29from os.path import dirname, basename
30from os import mkdir, access, F_OK
31from shutil import copy
32import time
33import sys
34
35# Related major packages
36from anuga.interface import create_domain_from_regions
37from anuga.interface import Dirichlet_boundary
38from anuga.interface import Reflective_boundary
39from anuga.interface import Field_boundary
40from anuga.interface import create_sts_boundary
41from anuga.interface import csv2building_polygons
42
43from anuga.shallow_water.data_manager import start_screen_catcher
44from anuga.shallow_water.data_manager import copy_code_files
45from anuga.utilities.polygon import read_polygon, Polygon_function
46   
47# Application specific imports
48import project  # Definition of file names and polygons
49
50   
51#-----------------------------------------------------------------------
52# Copy scripts to time stamped output directory and capture screen
53# output to file
54#-----------------------------------------------------------------------
55
56#copy script must be before screen_catcher
57#copy_code_files(project.output_run_time_dir, __file__,
58#         dirname(project.__file__)+sep+ project.__name__+'.py' )
59#start_screen_catcher(project.output_run_time_dir, myid, numprocs)
60
61
62#--------------------------------------------------------------------------
63# Create the computational domain based on overall clipping polygon with
64# a tagged boundary and interior regions defined in project.py along with
65# resolutions (maximal area of per triangle) for each polygon
66#--------------------------------------------------------------------------
67print 'Create computational domain'
68
69# Read in boundary from ordered sts file
70urs_bounding_polygon = create_sts_boundary(project.urs_boundary_name)
71
72# Reading the landward defined points, this incorporates the original clipping
73# polygon minus the 100m contour
74landward_bounding_polygon = read_polygon(project.landward_dir)
75
76# Combine sts polyline with landward points
77bounding_polygon = urs_bounding_polygon + landward_bounding_polygon
78
79# Number of boundary segments
80N = len(urs_bounding_polygon)-1
81
82# Boundary tags refer to project.landward 4 points equals 5 segments start at N
83boundary_tags={'back': [N+1,N+2,N+3,N+4, N+5],
84               'side': [N,N+6],
85               'ocean': range(N)}
86
87# Build mesh and domain
88domain = create_domain_from_regions(bounding_polygon,
89                                    boundary_tags=boundary_tags,
90                                    maximum_triangle_area=project.res_poly_all,
91                                    interior_regions=project.interior_regions,
92                                    mesh_filename=project.meshes_dir_name,
93                                    use_cache=True,
94                                    verbose=True)
95print domain.statistics()
96
97
98domain.set_name(project.scenario_name)
99domain.set_datadir(project.output_run_time_dir) 
100domain.set_minimum_storable_height(0.01)    # Don't store depth less than 1cm
101
102#-------------------------------------------------------------------------
103# Setup initial conditions
104#-------------------------------------------------------------------------
105print 'Setup initial conditions'
106
107# Set the initial stage in the offcoast region only
108IC = Polygon_function([(project.poly_mainland, 0),
109                       (project.poly_marina, 0)], 
110                      default=project.tide,
111                      geo_reference=domain.geo_reference)
112domain.set_quantity('stage', IC, use_cache=True, verbose=True)
113
114domain.set_quantity('friction', project.friction) 
115
116domain.set_quantity('elevation', 
117                    filename=project.combined_dir_name+'.pts',
118                    use_cache=True,
119                    verbose=True,
120                    alpha = project.alpha)
121
122
123#-------------------------------------------------------------------------
124# Setup boundary conditions
125#-------------------------------------------------------------------------
126
127print 'Set boundary - available tags:', domain.get_boundary_tags()
128
129boundary_urs_out=project.boundaries_dir_event + sep + project.scenario_name
130
131Br = Reflective_boundary(domain)
132Bd = Dirichlet_boundary([project.tide,0,0])
133
134
135Bf = Field_boundary(boundary_urs_out+'.sts',
136                    domain, mean_stage=project.tide,
137                    time_thinning=1,
138                    default_boundary=Bd,
139                    boundary_polygon=bounding_polygon,                   
140                    use_cache=True,
141                    verbose=True)
142
143
144domain.set_boundary({'back': Br,
145                     'side': Bd,
146                     'ocean': Bf}) 
147
148#-------------------------------------------------------------------------
149# Evolve system through time
150#-------------------------------------------------------------------------
151t0 = time.time()
152for t in domain.evolve(yieldstep=project.yieldstep, 
153                       finaltime=project.finaltime,
154                       skip_initial_step=False): 
155    print domain.timestepping_statistics()
156    print domain.boundary_statistics(tags='ocean')
157
158print 'Simulation took %.2f seconds' %(time.time()-t0)
159
160   
161   
162
163     
164   
Note: See TracBrowser for help on using the repository browser.