source: anuga_work/production/karratha_2006/run_karratha.py @ 3716

Last change on this file since 3716 was 3716, checked in by ole, 17 years ago

Ongoing work on run_okushiri and karratha_2006

File size: 4.7 KB
Line 
1"""Script for running tsunami inundation scenario for Karratha, 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
49
50
51# creates copy of code in output dir if dir doesn't exist
52if access(project.outputtimedir,F_OK) == 0:
53    mkdir (project.outputtimedir)
54copy (dirname(project.__file__) +sep+ project.__name__+'.py',
55      project.outputtimedir + project.__name__+'.py')
56copy (__file__, project.outputtimedir + basename(__file__))
57print 'project.outputtimedir',project.outputtimedir
58
59
60
61#--------------------------------------------------------------------------
62# Create the triangular mesh based on overall clipping polygon with a
63# tagged
64# boundary and interior regions defined in project.py along with
65# resolutions (maximal area of per triangle) for each polygon
66#--------------------------------------------------------------------------
67
68
69interior_regions = [#[project.karratha_polygon, 25000],
70                    #[project.dampier_polygon, 2000],
71                    #[project.refinery_polygon, 2000],
72                    #[project.point_polygon, 2000]]
73                    [project.neil1_polygon, 4000],
74                    [project.neil2_polygon, 64000]]
75
76
77
78print 'start create mesh from regions'
79meshname = project.meshname + '.msh'
80create_mesh_from_regions(project.bounding_polygon,
81                         boundary_tags={'back': [7, 8], 'side': [0, 6],
82                                        'ocean': [1, 2, 3, 4, 5]},
83                         maximum_triangle_area=200000,
84                         filename=meshname,
85                         use_cache=True,
86                         verbose=True)
87
88#-------------------------------------------------------------------------
89# Setup computational domain
90#-------------------------------------------------------------------------
91domain = Domain(meshname, use_cache=True, verbose=True)
92print domain.statistics()
93domain.set_name(project.basename)
94domain.set_datadir(project.outputtimedir)
95
96
97#-------------------------------------------------------------------------
98# Setup initial conditions
99#-------------------------------------------------------------------------
100tide = 0.
101domain.set_quantity('stage', tide)
102domain.set_quantity('friction', 0.0) 
103domain.set_quantity('elevation', 
104                    filename = project.datadir + project.basename + '.pts',
105                    use_cache = False,
106                    verbose = True,
107                    alpha = 0.1)
108
109#-------------------------------------------------------------------------
110# Setup boundary conditions
111#-------------------------------------------------------------------------
112
113print 'Available boundary tags', domain.get_boundary_tags()
114
115Bf = File_boundary(source_dir + project.boundary_basename + '.sww',
116                   domain, verbose = True)
117Br = Reflective_boundary(domain)
118Bd = Dirichlet_boundary([tide,0,0])
119domain.set_boundary({'back': Br,
120                     'side': Bd,
121                     'ocean': Bf}) 
122
123
124#----------------------------------------------------------------------------
125# Evolve system through time
126#----------------------------------------------------------------------------
127import time
128t0 = time.time()
129
130for t in domain.evolve(yieldstep = 60, finaltime = 28600):
131    domain.write_time()
132    domain.write_boundary_statistics(tags = 'ocean')     
133   
134print 'That took %.2f seconds' %(time.time()-t0)
135
Note: See TracBrowser for help on using the repository browser.