source: anuga_work/production/dampier_2006/run_dampier.py @ 3871

Last change on this file since 3871 was 3871, checked in by nick, 17 years ago

dampier updates

File size: 6.6 KB
Line 
1"""Script for running tsunami inundation scenario for Dampier, WA, 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 project.output_run_time_dir
6
7The scenario is defined by a triangular mesh created from project.polygon,
8the elevation data and a simulated submarine landslide.
9
10Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton, Nick Bartzis, GA - 2006
11"""
12
13#------------------------------------------------------------------------------
14# Import necessary modules
15#------------------------------------------------------------------------------
16
17# Standard modules
18from os import sep
19from os.path import dirname, basename
20from os import mkdir, access, F_OK
21from shutil import copy
22import time
23import sys
24
25
26# Related major packages
27from anuga.shallow_water import Domain
28from anuga.shallow_water import Dirichlet_boundary
29from anuga.shallow_water import File_boundary
30from anuga.shallow_water import Reflective_boundary
31
32from anuga.pmesh.mesh_interface import create_mesh_from_regions
33
34from anuga.geospatial_data.geospatial_data import *
35from anuga.abstract_2d_finite_volumes.util import Screen_Catcher
36
37# Application specific imports
38import project                 # Definition of file names and polygons
39
40
41
42#------------------------------------------------------------------------------
43# Copy scripts to time stamped output directory and capture screen
44# output to file
45#------------------------------------------------------------------------------
46
47# filenames
48
49build_time = '20061025_113524_build'
50boundaries_name = project.boundaries_name
51meshes_dir_name = project.meshes_dir_name+'.msh'
52#source_dir = project.boundarydir
53boundaries_time_dir_name = project.boundaries_dir + build_time + sep + boundaries_name
54
55#from anuga.shallow_water.data_manager import urs2sww
56
57# creates copy of code in output dir if dir doesn't exist
58if access(project.output_run_time_dir,F_OK) == 0:
59    print 'project.output_run_time_dir',dirname(project.output_run_time_dir)
60    mkdir (project.output_run_time_dir,0777)
61copy (dirname(project.__file__) +sep+ project.__name__+'.py',
62      project.output_run_time_dir + project.__name__+'.py')
63copy (__file__, project.output_run_time_dir + basename(__file__))
64print 'project.output_run_time_dir',project.output_run_time_dir
65
66#normal screen output is stored in
67screen_output_name = project.output_run_time_dir + "screen_output.txt"
68screen_error_name = project.output_run_time_dir + "screen_error.txt"
69
70#used to catch screen output to file
71sys.stdout = Screen_Catcher(screen_output_name)
72sys.stderr = Screen_Catcher(screen_error_name)
73
74print 'USER:    ', project.user
75
76#--------------------------------------------------------------------------
77# Create the triangular mesh based on overall clipping polygon with a
78# tagged
79# boundary and interior regions defined in project.py along with
80# resolutions (maximal area of per triangle) for each polygon
81#--------------------------------------------------------------------------
82
83interior_regions = [#[project.karratha_polygon, 25000],
84                    #[project.dampier_polygon, 2000],
85                    #[project.refinery_polygon, 2000],
86                    #[project.point_polygon, 2000]]
87                    #[project.cipma_polygon, 20000]]
88                    [project.cipma_polygon, 50000]]    # Caused memory error?
89
90#---------------------------
91# this is check that no interior polygon is outside the bounding poly
92#------------------------------
93
94"""
95count = 0
96for i in range(len(interior_regions)):
97    region = interior_regions[i]
98    interior_polygon = region[0]
99    if len(inside_polygon(interior_polygon, bounding_polygon,
100                   closed = True, verbose = False)) <> len(interior_polygon):
101        print 'WARNING: interior polygon %d is outside bounding polygon' %(i)
102        count += 1
103
104if count == 0:
105    print 'interior regions OK'
106else:
107    print 'check out your interior polygons'
108    print 'check %s in production directory' %figname
109    import sys; sys.exit()
110"""
111#-------------------------------------
112
113print 'start create mesh from regions'
114meshes_dir_name = project.meshes_dir_name + '.msh'
115create_mesh_from_regions(project.bounding_polygon,
116                         boundary_tags={'back': [7, 8], 'side': [0, 6],
117                                        'ocean': [1, 2, 3, 4, 5]},
118                         maximum_triangle_area=300000,
119                         interior_regions=interior_regions,
120                         filename=meshes_dir_name,
121                         use_cache=True,
122                         verbose=True)
123
124#-------------------------------------------------------------------------
125# Setup computational domain
126#-------------------------------------------------------------------------
127print 'Setup computational domain'
128domain = Domain(meshes_dir_name, use_cache=True, verbose=True)
129print domain.statistics()
130domain.set_name(project.scenario_name)
131domain.set_datadir(project.output_run_time_dir)
132domain.set_default_order(2)
133domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
134
135
136
137#-------------------------------------------------------------------------
138# Setup initial conditions
139#-------------------------------------------------------------------------
140tide = project.tide
141domain.set_quantity('stage', tide)
142domain.set_quantity('friction', 0.0) 
143#combined_time_dir_name = project.topographies_dir+build_time+project.combined_name
144domain.set_quantity('elevation', 
145                    filename = project.topographies_dir + build_time + sep + project.combined_name + '.pts',
146                    use_cache = True,
147                    verbose = True,
148                    alpha = 0.1)
149
150
151#-------------------------------------------------------------------------
152# Setup boundary conditions
153#-------------------------------------------------------------------------
154
155
156print 'Available boundary tags', domain.get_boundary_tags()
157
158
159Bf = File_boundary(boundaries_time_dir_name + '.sww',
160                   domain, verbose = True)
161#Bf = File_boundary(source_dir + project.boundary_scenario_name + '.sww',
162#                   domain, verbose = True)
163Br = Reflective_boundary(domain)
164Bd = Dirichlet_boundary([tide,0,0])
165domain.set_boundary({'back': Br,
166                     'side': Bd,
167                     'ocean': Bf}) 
168
169#----------------------------------------------------------------------------
170# Evolve system through time
171#----------------------------------------------------------------------------
172import time
173t0 = time.time()
174
175for t in domain.evolve(yieldstep = 120, finaltime = 28800):
176    domain.write_time()
177    domain.write_boundary_statistics(tags = 'ocean')     
178   
179print 'That took %.2f seconds' %(time.time()-t0)
180
Note: See TracBrowser for help on using the repository browser.