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

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

build_dampier.py is working
run_dampier.py should work... it is still running
project.py is working

File size: 6.2 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 *
35
36# Application specific imports
37import project                 # Definition of file names and polygons
38
39
40
41#------------------------------------------------------------------------------
42# Copy scripts to time stamped output directory and capture screen
43# output to file
44#------------------------------------------------------------------------------
45
46# filenames
47
48build_time = '20061025_113524_build'
49boundaries_name = project.boundaries_name
50meshes_dir_name = project.meshes_dir_name+'.msh'
51#source_dir = project.boundarydir
52boundaries_time_dir_name = project.boundaries_dir + build_time + sep + boundaries_name
53
54#from anuga.shallow_water.data_manager import urs2sww
55
56# creates copy of code in output dir if dir doesn't exist
57if access(project.output_run_time_dir,F_OK) == 0:
58    print 'project.output_run_time_dir',dirname(project.output_run_time_dir)
59    mkdir (project.output_run_time_dir,0777)
60copy (dirname(project.__file__) +sep+ project.__name__+'.py',
61      project.output_run_time_dir + project.__name__+'.py')
62copy (__file__, project.output_run_time_dir + basename(__file__))
63print 'project.output_run_time_dir',project.output_run_time_dir
64
65#--------------------------------------------------------------------------
66# Create the triangular mesh based on overall clipping polygon with a
67# tagged
68# boundary and interior regions defined in project.py along with
69# resolutions (maximal area of per triangle) for each polygon
70#--------------------------------------------------------------------------
71
72interior_regions = [#[project.karratha_polygon, 25000],
73                    #[project.dampier_polygon, 2000],
74                    #[project.refinery_polygon, 2000],
75                    #[project.point_polygon, 2000]]
76                    #[project.cipma_polygon, 20000]]
77                    [project.cipma_polygon, 50000]]    # Caused memory error?
78
79#---------------------------
80# this is check that no interior polygon is outside the bounding poly
81#------------------------------
82
83"""
84count = 0
85for i in range(len(interior_regions)):
86    region = interior_regions[i]
87    interior_polygon = region[0]
88    if len(inside_polygon(interior_polygon, bounding_polygon,
89                   closed = True, verbose = False)) <> len(interior_polygon):
90        print 'WARNING: interior polygon %d is outside bounding polygon' %(i)
91        count += 1
92
93if count == 0:
94    print 'interior regions OK'
95else:
96    print 'check out your interior polygons'
97    print 'check %s in production directory' %figname
98    import sys; sys.exit()
99"""
100#-------------------------------------
101
102print 'start create mesh from regions'
103meshes_dir_name = project.meshes_dir_name + '.msh'
104create_mesh_from_regions(project.bounding_polygon,
105                         boundary_tags={'back': [7, 8], 'side': [0, 6],
106                                        'ocean': [1, 2, 3, 4, 5]},
107                         maximum_triangle_area=300000,
108                         interior_regions=interior_regions,
109                         filename=meshes_dir_name,
110                         use_cache=True,
111                         verbose=True)
112
113#-------------------------------------------------------------------------
114# Setup computational domain
115#-------------------------------------------------------------------------
116print 'Setup computational domain'
117domain = Domain(meshes_dir_name, use_cache=True, verbose=True)
118print domain.statistics()
119domain.set_name(project.scenario_name)
120domain.set_datadir(project.output_run_time_dir)
121domain.set_default_order(2)
122domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
123
124
125
126#-------------------------------------------------------------------------
127# Setup initial conditions
128#-------------------------------------------------------------------------
129tide = 2.4
130domain.set_quantity('stage', tide)
131domain.set_quantity('friction', 0.0) 
132#combined_time_dir_name = project.topographies_dir+build_time+project.combined_name
133domain.set_quantity('elevation', 
134                    filename = project.topographies_dir + build_time + sep + project.combined_name + '.pts',
135                    use_cache = True,
136                    verbose = True,
137                    alpha = 0.1)
138
139
140#-------------------------------------------------------------------------
141# Setup boundary conditions
142#-------------------------------------------------------------------------
143
144
145print 'Available boundary tags', domain.get_boundary_tags()
146
147
148Bf = File_boundary(boundaries_time_dir_name + '.sww',
149                   domain, verbose = True)
150#Bf = File_boundary(source_dir + project.boundary_scenario_name + '.sww',
151#                   domain, verbose = True)
152Br = Reflective_boundary(domain)
153Bd = Dirichlet_boundary([tide,0,0])
154domain.set_boundary({'back': Br,
155                     'side': Bd,
156                     'ocean': Bf}) 
157
158#----------------------------------------------------------------------------
159# Evolve system through time
160#----------------------------------------------------------------------------
161import time
162t0 = time.time()
163
164for t in domain.evolve(yieldstep = 120, finaltime = 28800):
165    domain.write_time()
166    domain.write_boundary_statistics(tags = 'ocean')     
167   
168print 'That took %.2f seconds' %(time.time()-t0)
169
Note: See TracBrowser for help on using the repository browser.