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

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

updates to dampier

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