source: anuga_work/production/australia_ph2/ceduna/run_model.py @ 6398

Last change on this file since 6398 was 6398, checked in by myall, 16 years ago

ceduna working using older scripts - going to upgrade to the latest ones anyway, so this isn't very important

File size: 5.9 KB
Line 
1"""Run a tsunami inundation scenario for Busselton, WA, Australia.
2
3The scenario is defined by a triangular mesh created from project.polygon, the
4elevation data is compiled into a pts file through build_elevation.py and a
5simulated 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_elevation.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#------------------------------------------------------------------------------
19# Import necessary modules
20#------------------------------------------------------------------------------
21
22# Standard modules
23import os
24import time
25
26# Related major packages
27from anuga.interface import create_domain_from_regions
28from anuga.interface import Dirichlet_boundary
29from anuga.interface import Reflective_boundary
30from anuga.interface import Field_boundary
31from anuga.interface import create_sts_boundary
32from anuga.interface import csv2building_polygons
33from file_length import file_length
34
35from anuga.shallow_water.data_manager import start_screen_catcher
36from anuga.shallow_water.data_manager import copy_code_files
37from anuga.utilities.polygon import read_polygon, Polygon_function
38   
39# Application specific imports
40import project  # Definition of file names and polygons
41
42
43#------------------------------------------------------------------------------
44# Copy scripts to time stamped output directory and capture screen
45# output to file. Copy script must be before screen_catcher
46#------------------------------------------------------------------------------
47copy_code_files(project.output_run, __file__, 
48                os.path.dirname(project.__file__)+os.sep+\
49                project.__name__+'.py' )
50start_screen_catcher(project.output_run, 0, 1)
51
52
53#------------------------------------------------------------------------------
54# Create the computational domain based on overall clipping polygon with
55# a tagged boundary and interior regions defined in project.py along with
56# resolutions (maximal area of per triangle) for each polygon
57#------------------------------------------------------------------------------
58print 'Create computational domain'
59
60# Read in boundary from ordered sts file
61event_sts = create_sts_boundary(project.event_sts)
62print 'HELLO event_sts', event_sts
63
64# Reading the landward defined points, this incorporates the original clipping
65# polygon minus the 100m contour
66landward_boundary = read_polygon(project.landward_boundary)
67print 'landward boundary', landward_boundary
68
69# Combine sts polyline with landward points
70bounding_polygon_sts = event_sts + landward_boundary
71print 'bounding polygon', bounding_polygon_sts
72
73# Number of boundary segments
74N = len(event_sts)-1
75# Number of landward_boundary points
76M = file_length(project.landward_boundary)
77
78# Boundary tags refer to project.landward_boundary
79# 4 points equals 5 segments start at N
80boundary_tags={'back': range(N+1,N+M),
81               'side': [N,N+M],
82               'ocean': range(N)}
83
84# Build mesh and domain
85domain = create_domain_from_regions(bounding_polygon_sts,
86                                    boundary_tags=boundary_tags,
87                                    maximum_triangle_area=project.bounding_maxarea,
88                                    interior_regions=project.interior_regions,
89                                    mesh_filename=project.meshes,
90                                    use_cache=False,
91                                    verbose=False)
92print domain.statistics()
93
94domain.set_name(project.scenario_name)
95domain.set_datadir(project.output_run) 
96domain.set_minimum_storable_height(0.01)    # Don't store depth less than 1cm
97
98
99#------------------------------------------------------------------------------
100# Setup initial conditions
101#------------------------------------------------------------------------------
102print 'Setup initial conditions'
103
104# Set the initial stage in the offcoast region only
105##IC = Polygon_function(project.land_initial_conditions,
106##                      default=project.tide,
107##                      geo_reference=domain.geo_reference)
108domain.set_quantity('stage', 0, use_cache=True, verbose=True)
109domain.set_quantity('friction', project.friction) 
110domain.set_quantity('elevation', 
111                    filename=project.combined_elevation+'.pts',
112                    use_cache=True,
113                    verbose=True,
114                    alpha=project.alpha)
115
116
117#------------------------------------------------------------------------------
118# Setup boundary conditions
119#------------------------------------------------------------------------------
120print 'Set boundary - available tags:', domain.get_boundary_tags()
121
122Br = Reflective_boundary(domain)
123Bd = Dirichlet_boundary([project.tide,0,0])
124Bf = Field_boundary(project.event_sts+'.sts',
125                    domain, mean_stage=project.tide,
126                    time_thinning=1,
127                    default_boundary=Bd,
128                    boundary_polygon=bounding_polygon_sts,                   
129                    use_cache=True,
130                    verbose=True)
131
132
133domain.set_boundary({'back': Br,
134                     'side': Bd,
135                     'ocean': Bf}) 
136
137
138#------------------------------------------------------------------------------
139# Evolve system through time
140#------------------------------------------------------------------------------
141t0 = time.time()
142for t in domain.evolve(yieldstep=project.yieldstep, 
143                       finaltime=project.finaltime,
144                       skip_initial_step=False): 
145    print domain.timestepping_statistics()
146    print domain.boundary_statistics(tags='ocean')
147
148print 'Simulation took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.