source: production/karratha_2005/run_karratha.py @ 1920

Last change on this file since 1920 was 1920, checked in by ole, 19 years ago

Karratha

File size: 5.0 KB
Line 
1"""Script for running a 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.outputdir
6
7The scenario is defined by a triangular mesh created from project.polygon,
8the elevation data and boundary data obtained from a tsunami simulation done with MOST.
9
10Ole Nielsen, GA - 2005
11"""
12
13#tide = 0.75   #HMWS estimate by Colin French, GA
14tide = 0       #HMW
15
16
17import os
18import time
19
20
21from pyvolution.shallow_water import Domain, Reflective_boundary, File_boundary,\
22     Dirichlet_boundary, Time_boundary, Transmissive_boundary
23from pyvolution.data_manager import convert_dem_from_ascii2netcdf,\
24     dem2pts, ferret2sww
25from pyvolution.pmesh2domain import pmesh_to_domain_instance
26from caching import cache
27import project
28
29#Data preparation
30#Convert ASC 2 DEM 2 PTS using source data and store result in source data
31demname = project.demname
32
33cache(convert_dem_from_ascii2netcdf, demname, {'verbose': True},
34      dependencies = [demname + '.asc'],
35      verbose = True)
36      #evaluate = True)
37
38cache(dem2pts, demname, {'verbose': True},
39      dependencies = [demname + '.dem'],     
40      verbose = True)
41
42
43#Convert MOST boundary
44source_dir = project.boundarydir
45
46from pyvolution.data_manager import ferret2sww
47
48south = project.south
49north = project.north
50west = project.west
51east = project.east
52
53
54cache(ferret2sww,
55      (source_dir+project.boundary_basename,
56       project.boundary_basename), 
57      {'verbose': True,
58       'minlat': south-1,
59       'maxlat': north+1,
60       'minlon': west-1,
61       'maxlon': east+1,
62       'origin': project.mesh_origin,
63       'mean_stage': tide,
64       'zscale': 1,                 #Enhance tsunami
65       'fail_on_NaN': False,
66       'inverted_bathymetry': True},
67      #evaluate = True,
68      verbose = True)
69#FIXME: Dependencies
70
71   
72#ferret2sww(source_dir+project.boundary_basename,
73#           project.boundary_basename,
74#           verbose=True,
75#           minlat=south-1, maxlat=north+1,
76#           minlon=west-1, maxlon=east+1,
77#           origin = project.mesh_origin,
78#           mean_stage = tide,
79#           zscale = 1,
80#           fail_on_NaN = False,
81#           inverted_bathymetry = True)
82
83
84#Create Triangular Mesh
85from create_mesh import create_mesh
86
87interior_regions = [[project.karratha_polygon, 25000],
88                    [project.dampier_polygon, 8000],
89                    [project.refinery_polygon, 8000]]
90
91m = cache(create_mesh,
92          project.polygon,
93          {'boundary_tags': {'back': [7, 8], 'side': [0, 6], 'ocean': [1, 2, 3, 4, 5]},
94          'resolution': 100000,
95          'filename': project.meshname + '.msh',
96          'interior_regions': interior_regions},
97          verbose = True)         
98          #verbose = True,
99          #evaluate = True )
100
101
102#Setup domain
103mesh = project.meshname + '.msh'
104domain = cache(pmesh_to_domain_instance, (mesh, Domain),
105               dependencies = [mesh],                     
106               verbose = True)               
107
108
109domain.set_name(project.basename)
110domain.set_datadir(project.outputdir)
111domain.store = True
112
113domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
114#domain.quantities_to_be_stored = ['stage']
115       
116print 'Number of triangles = ', len(domain)
117print 'The extent is ', domain.get_extent()
118
119
120
121#Setup Initial Conditions
122domain.set_quantity('friction', 0)
123domain.set_quantity('stage', tide)
124domain.set_quantity('elevation',
125                    filename = demname + '.pts',
126                    use_cache = True,
127                    verbose = True)
128
129
130
131#Setup Boundary Conditions
132print domain.get_boundary_tags()
133
134Bf = File_boundary(project.boundary_basename + '.sww', domain, verbose = True)
135#domain.starttime = 3000  #Obtained from MOST
136domain.starttime = 0  #Obtained from MOST
137
138Br = Reflective_boundary(domain)
139Bt = Transmissive_boundary(domain)
140Bd = Dirichlet_boundary([tide,0,0])
141Bw = Time_boundary(domain=domain,
142                   f=lambda t: [(60<t<660)*4, 0, 0])
143
144
145domain.set_boundary( {'back': Br,'side': Bd, 'ocean': Bf} ) #MOST tsunami
146
147#domain.set_boundary( {'back': Br,'side': Bd, 'ocean': Bd} )  #Nothing
148
149
150#Evolve
151import time
152t0 = time.time()
153
154#for t in domain.evolve(yieldstep = 60, finaltime = 15000):
155#    domain.write_time()
156#    domain.write_boundary_statistics(tags = 'ocean') #quantities = 'stage')
157
158#for t in domain.evolve(yieldstep = 20, finaltime = 35000,
159#                       skip_initial_step = True):
160#    domain.write_time()
161#    domain.write_boundary_statistics(tags = 'ocean') #quantities = 'stage')
162
163for t in domain.evolve(yieldstep = 60, finaltime = 40000):
164                       #skip_initial_step = True):
165    domain.write_time()
166    domain.write_boundary_statistics(tags = 'ocean') #quantities = 'stage')       
167
168print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.