source: production/onslow_2006/run_onslow_old.py @ 3514

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

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

File size: 9.7 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#------------------------------------------------------------------------------
15# Import necessary modules
16#------------------------------------------------------------------------------
17
18# Standard modules
19import os
20import time
21
22# Related major packages
23from anuga.pyvolution.shallow_water import Domain, Reflective_boundary, \
24                            Dirichlet_boundary, Time_boundary, File_boundary
25from anuga.pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts
26from anuga.pyvolution.combine_pts import combine_rectangular_points_files
27from anuga.pyvolution.pmesh2domain import pmesh_to_domain_instance
28#from anuga.geospatial_data.geospatial_data import *
29
30# Application specific imports
31import project                 # Definition of file names and polygons
32from smf import slump_tsunami  # Function for submarine mudslide
33
34from shutil import copy
35from os import mkdir, access, F_OK
36
37#------------------------------------------------------------------------------
38# Preparation of topographic data
39#
40# Convert ASC 2 DEM 2 PTS using source data and store result in source data
41# Do for coarse and fine data
42# Fine pts file to be clipped to area of interest
43#------------------------------------------------------------------------------
44
45# filenames
46coarsedemname = project.coarsedemname
47
48onshore_dem_name = project.onshore_dem_name
49
50meshname = project.meshname+'.msh'
51
52#############################################################
53# making a copy of project and run to a time stamped directory
54# with the ouput
55source_dir = project.boundarydir
56
57#if dir doesn't exists then makes dir
58if access(project.outputdir,F_OK) == 0 :
59    mkdir (project.outputdir)
60# creates copy of code in output dir
61copy (project.codedirname, project.outputdir + project.codename)
62copy (project.codedir + 'run_onslow.py', project.outputdir + 'run_onlsow.py')
63#print "copied to"+ project.outputdir + project.codename + 'and run_onlsow.py'
64#################################################################
65'''
66# coarse data
67convert_dem_from_ascii2netcdf(coarsedemname, use_cache=True, verbose=True)
68dem2pts(coarsedemname, use_cache=True, verbose=True)
69
70
71# fine data (clipping the points file to smaller area)
72convert_dem_from_ascii2netcdf(onshore_dem_name, use_cache=True, verbose=True)
73dem2pts(onshore_dem_name,
74        easting_min=project.eastingmin,
75        easting_max=project.eastingmax,
76        northing_min=project.northingmin,
77        northing_max= project.northingmax,
78        use_cache=True,
79        verbose=True)
80
81print 'before off xya to object'
82offshore_pts = Geospatial_data(project.offshore_dem_name + '.xya')
83print 'before offshore to dict'
84offshore_dict = geospatial_data2points_dictionary(offshore_pts)
85
86print 'offshore to pts file'
87export_points_file(project.offshore_dem_name + '.pts', offshore_dict)
88'''
89# combining the coarse and fine data
90# NOTE MUST HAVE FINE FIRST!
91'''
92combine_rectangular_points_files(
93                                 project.onshore_dem_name + '.pts',
94                                 project.coarsedemname + '.pts',
95                                 project.combined_dem_name + '.pts')
96
97print 'create G1'
98G1 = Geospatial_data(project.onshore_dem_name + '.pts')
99print 'G1 dict'
100G1_points_dict = geospatial_data2points_dictionary(G1)
101print 'G1 xya file'
102export_points_file(project.datadir + 'offshore_dem.xya', G1_points_dict)
103print 'create G2'
104G2 = Geospatial_data(project.offshore_dem_name + '.xya')
105print 'G1+g2'
106G = G1 + G2
107print 'G dict'
108G_points_dict = geospatial_data2points_dictionary(G)
109print 'G to xya'
110export_points_file(project.combined_dem_name + '.xya', G_points_dict)
111
112
113add_points_files(project.onshore_dem_name + '.pts',
114                  project.offshore_dem_name + '.xya',
115                  project.combined_dem_name + '.pts')
116
117print 'please'
118
119add_points_files(project.onshore_dem_name + '.pts',
120                  project.offshore_dem_name + '.pts',
121                  project.combined_dem_name + '.pts')
122
123print ' finished with points'
124'''
125#-------------------------------------------------------------------------------                                 
126# Create the triangular mesh based on overall clipping polygon with a tagged
127# boundary and interior regions defined in project.py along with
128# resolutions (maximal area of per triangle) for each polygon
129#------------------------------------------------------------------------------
130from pmesh.mesh_interface import create_mesh_from_regions
131
132# original
133interior_res = 50000
134interior_regions = [[project.poly_onslow, interior_res],
135                    [project.poly_thevenard, interior_res],
136                    [project.poly_direction, interior_res]]
137                    #[project.testpoly, interior_res]]
138print 'number of interior regions', len(interior_regions)
139
140from caching import cache
141_ = cache(create_mesh_from_regions,
142          project.polyAll,
143          {'boundary_tags': {'top': [0], 'topleft': [1],
144                             'left': [2], 'bottom': [3],
145                             'bottomright': [4], 'topright': [5]},
146           'maximum_triangle_area': 1000000,
147           'filename': meshname,           
148           'interior_regions': interior_regions},
149          verbose = True)
150
151
152#------------------------------------------------------------------------------
153# Setup computational domain
154#------------------------------------------------------------------------------
155
156domain = pmesh_to_domain_instance(meshname, Domain,
157                                  use_cache = True,
158                                  verbose = True)
159
160print 'Number of triangles = ', len(domain)
161print 'The extent is ', domain.get_extent()
162print domain.statistics()
163
164domain.set_name(project.basename)
165domain.set_datadir(project.outputdir)
166domain.set_quantities_to_be_stored(['stage'])
167
168print 'hi'
169#------------------------------------------------------------------------------
170# Set up scenario (tsunami_source is a callable object used with set_quantity)
171#------------------------------------------------------------------------------
172'''
173tsunami_source = slump_tsunami(length=30000.0,
174                               depth=400.0,
175                               slope=6.0,
176                               thickness=176.0,
177                               radius=3330,
178                               dphi=0.23,
179                               x0=project.slump_origin[0],
180                               y0=project.slump_origin[1],
181                               alpha=0.0,
182                               domain=domain)
183
184'''
185#------------------------------------------------------------------------------
186# Setup initial conditions
187#------------------------------------------------------------------------------
188
189tide = 0.
190
191domain.set_quantity('stage', tide)
192domain.set_quantity('friction', 0.0) 
193print 'hi1', project.combined_dem_name + '.pts'
194domain.set_quantity('elevation', 
195#                    0.
196#                    filename = project.onshore_dem_name + '.pts',
197                    filename = project.combined_dem_name + '.pts',
198#                    filename = project.coarsedemname + '.pts',
199                    use_cache = True,
200                    verbose = True
201                    )
202
203print 'hi2'
204#------------------------------------------------------------------------------
205# Setup boundary conditions (all reflective)
206#------------------------------------------------------------------------------
207'''
208from anuga.pyvolution.data_manager import ferret2sww
209
210south = project.south
211north = project.north
212west = project.west
213east = project.east
214
215cache(ferret2sww,
216      (source_dir + project.boundary_basename,
217       source_dir + project.boundary_basename),
218      {'verbose': True,
219# note didn't work with the below
220#       'minlat': south - 1,
221#       'maxlat': north + 1,
222#       'minlon': west - 1,
223#       'maxlon': east + 1,
224       'minlat': south,
225       'maxlat': north,
226       'minlon': west,
227       'maxlon': east,
228#       'origin': project.mesh_origin,
229       'origin': domain.geo_reference.get_origin(),
230       'mean_stage': tide,
231       'zscale': 1,                 #Enhance tsunami
232       'fail_on_NaN': False,
233       'inverted_bathymetry': True},
234      #evaluate = True,
235       verbose = True)
236'''
237
238print 'Available boundary tags', domain.get_boundary_tags()
239
240#Bf = File_boundary(source_dir + project.boundary_basename + '.sww',
241#                    domain, verbose = True)
242Br = Reflective_boundary(domain)
243Bd = Dirichlet_boundary([tide,0,0])
244
245print 'hi3'
246# 7 min square wave starting at 1 min, 6m high
247Bw = Time_boundary(domain = domain,
248                   f=lambda t: [(60<t<480)*6, 0, 0])
249
250domain.set_boundary( {'top': Bw, 'topleft': Bw,
251                             'left': Br, 'bottom': Br,
252                             'bottomright': Br, 'topright': Br} )
253print 'hi4'
254
255#------------------------------------------------------------------------------
256# Evolve system through time
257#------------------------------------------------------------------------------
258
259import time
260t0 = time.time()
261
262for t in domain.evolve(yieldstep = 50, finaltime = 50): 
263    domain.write_time()
264    domain.write_boundary_statistics(tags = 'top')     
265   
266print 'That took %.2f seconds' %(time.time()-t0)
267print 'hi5'
Note: See TracBrowser for help on using the repository browser.