source: production/MOST_example/run_onslow.py @ 3553

Last change on this file since 3553 was 3553, checked in by duncan, 18 years ago

bug fixes

File size: 8.0 KB
Line 
1"""Script for running a tsunami inundation scenario for Onslow, 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 a simulated submarine landslide.
9
10Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis, GA - 2006
11"""
12
13
14#-------------------------------------------------------------------------------# Import necessary modules
15#-------------------------------------------------------------------------------
16
17# Standard modules
18import os
19import time
20
21# Related major packages
22from anuga.pyvolution.shallow_water import Domain, Reflective_boundary, \
23                            Dirichlet_boundary, Time_boundary, File_boundary
24from anuga.pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts
25from anuga.pyvolution.combine_pts import combine_rectangular_points_files
26from anuga.pyvolution.pmesh2domain import pmesh_to_domain_instance
27from anuga.fit_interpolate.fit import fit_to_mesh_file
28
29# Application specific imports
30import project                 # Definition of file names and polygons
31from smf import slump_tsunami  # Function for submarine mudslide
32
33from shutil import copy
34from os import mkdir, access, F_OK
35
36from anuga.geospatial_data.geospatial_data import *
37
38#-------------------------------------------------------------------------------
39# Preparation of topographic data
40#
41# Convert ASC 2 DEM 2 PTS using source data and store result in source data
42# Do for coarse and fine data
43# Fine pts file to be clipped to area of interest
44#-------------------------------------------------------------------------------
45
46# filenames
47coarsedemname = project.coarsedemname
48
49onshore_dem_name = project.onshore_dem_name
50
51offshore_points = project.offshore_dem_name
52
53meshname = project.meshname+'.msh'
54
55source_dir = project.boundarydir
56
57# creates copy of code in output dir
58if access(project.outputdir,F_OK) == 0 :
59    mkdir (project.outputdir)
60copy (project.codedirname, project.outputdir + project.codename)
61copy (project.codedir + 'run_onslow.py', project.outputdir + 'run_onslow.py')
62
63
64'''
65# coarse data
66convert_dem_from_ascii2netcdf(coarsedemname, use_cache=True, verbose=True)
67dem2pts(coarsedemname, use_cache=True, verbose=True)
68
69
70# fine data (clipping the points file to smaller area)
71convert_dem_from_ascii2netcdf(onshore_dem_name, use_cache=True, verbose=True)
72dem2pts(onshore_dem_name,
73        easting_min=project.eastingmin,
74        easting_max=project.eastingmax,
75        northing_min=project.northingmin,
76        northing_max= project.northingmax,
77        use_cache=True,
78        verbose=True)
79
80        '''
81print'create G1'
82G1 = Geospatial_data(file_name = project.offshore_dem_name + '.xya')
83
84print'create G2'
85G2 = Geospatial_data(file_name = project.onshore_dem_name + '.pts')
86
87print'add G1+G2'
88G = G1 + G2
89
90print'export G'
91G.new_export_points_file(project.combined_dem_name + '.pts')
92
93
94#-------------------------------------------------------------------------------                                 
95# Create the triangular mesh based on overall clipping polygon with a tagged
96# boundary and interior regions defined in project.py along with
97# resolutions (maximal area of per triangle) for each polygon
98#-------------------------------------------------------------------------------
99
100from anuga.pmesh.mesh_interface import create_mesh_from_regions
101
102# original
103interior_res = 1000000
104interior_regions = [[project.poly_onslow, interior_res],
105                    [project.poly_thevenard, interior_res],
106                    [project.poly_direction, interior_res]]
107                    #[project.testpoly, interior_res]]
108print 'number of interior regions', len(interior_regions)
109
110from caching import cache
111_ = cache(create_mesh_from_regions,
112          project.polyAll,
113          {'boundary_tags': {'top': [0], 'topleft': [1],
114                             'left': [2], 'bottom': [3],
115                             'bottomright': [4], 'topright': [5]},
116           'maximum_triangle_area': 10000000,
117           'filename': meshname,           
118           'interior_regions': interior_regions},
119          verbose = True)
120
121mesh_elevname = meshname[:-4] + '_elv' + meshname[-4:]
122
123cache(fit_to_mesh_file,(meshname,
124                 project.combined_dem_name + '.pts',
125                 mesh_elevname,
126                 DEFAULT_ALPHA),
127      {'verbose': True}
128      ,dependencies = [meshname, pointname]
129      #,evaluate = True     
130      ,verbose = False
131      )
132
133#-------------------------------------------------------------------------------                                 
134# Setup computational domain
135#-------------------------------------------------------------------------------                                 
136
137domain = pmesh_to_domain_instance(mesh_elevname, Domain,
138                                  use_cache = True,
139                                  verbose = True)
140
141print 'Number of triangles = ', len(domain)
142print 'The extent is ', domain.get_extent()
143print domain.statistics()
144
145domain.set_name(project.basename)
146domain.set_datadir(project.outputdir)
147domain.set_quantities_to_be_stored(['stage'])
148
149
150#-------------------------------------------------------------------------------                                 
151# Setup initial conditions
152#-------------------------------------------------------------------------------
153
154tide = 0.
155
156domain.set_quantity('stage', tide)
157domain.set_quantity('friction', 0.0) 
158print 'hi and file',project.combined_dem_name + '.pts'
159#domain.set_quantity('elevation',
160#                    0.
161#                    filename = project.onshore_dem_name + '.pts',
162#                    filename = project.combined_dem_name + '.pts',
163#                    filename = project.offshore_dem_name + '.pts',
164 #                   use_cache = False,
165  #                  verbose = True,
166   #                 alpha = 0.1
167    #                )
168print 'hi1'
169
170#-------------------------------------------------------------------------------                                 
171# Setup boundary conditions (all reflective)
172#-------------------------------------------------------------------------------
173
174from anuga.pyvolution.data_manager import ferret2sww
175
176south = project.south
177north = project.north
178west = project.west
179east = project.east
180
181cache(ferret2sww,
182      (source_dir + project.boundary_basename,
183       source_dir + project.boundary_basename), 
184      {'verbose': True,
185       'minlat': south,
186       'maxlat': north,
187       'minlon': west,
188       'maxlon': east,
189       'origin': domain.geo_reference.get_origin(),
190       'mean_stage': tide,
191       'zscale': 1,                 #Enhance tsunami
192       'fail_on_NaN': False,
193       'inverted_bathymetry': True},
194      #evaluate = True,
195       verbose = True)
196
197
198print 'Available boundary tags', domain.get_boundary_tags()
199
200Bf = File_boundary(source_dir + project.boundary_basename + '.sww', 
201                    domain, verbose = True)
202Br = Reflective_boundary(domain)
203Bd = Dirichlet_boundary([tide,0,0])
204
205
206# 7 min square wave starting at 1 min, 6m high
207Bw = Time_boundary(domain = domain,
208                   f=lambda t: [(60<t<480)*6, 0, 0])
209
210domain.set_boundary( {'top': Bf, 'topleft': Bf,
211                             'left': Br, 'bottom': Br,
212                             'bottomright': Br, 'topright': Bf} )
213
214
215#-------------------------------------------------------------------------------                                 
216# Evolve system through time
217#-------------------------------------------------------------------------------
218import time
219t0 = time.time()
220
221for t in domain.evolve(yieldstep = 50, finaltime = 100): 
222    domain.write_time()
223    domain.write_boundary_statistics(tags = 'top')     
224   
225print 'That took %.2f seconds' %(time.time()-t0)
226
227print 'finished'
Note: See TracBrowser for help on using the repository browser.