source: anuga_work/development/gong_2008/run_gong_slide_full.py @ 7840

Last change on this file since 7840 was 5677, checked in by sexton, 16 years ago

scripts ready for Bridgette

File size: 9.3 KB
Line 
1"""Script for running a tsunami inundation scenario for Wollongong, NSW, Australia.
2
3Source data such as elevation and boundary data is assumed to be available in
4directories specified by project_slide.py
5The output sww file is stored in project_slide.outputtimedir
6
7The scenario is defined by a triangular mesh created from project.polygon,
8the elevation data and a tsunami wave generated by s submarine mass failure.
9
10Bridgette Lewis, 2008
11"""
12
13#-------------------------------------------------------------------------------
14# Import necessary modules
15#-------------------------------------------------------------------------------
16
17# Standard modules
18import os
19import time
20from shutil import copy
21from os.path import dirname, basename
22from os import mkdir, access, F_OK, sep
23import sys
24
25# Related major packages
26from anuga.shallow_water import Domain, Reflective_boundary, Dirichlet_boundary
27from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
28from anuga.geospatial_data.geospatial_data import *
29from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
30
31# Application specific imports
32import project_slide              # Definition of file names and polygons
33
34#-------------------------------------------------------------------------------
35# Copy scripts to time stamped output directory and capture screen
36# output to file
37#-------------------------------------------------------------------------------
38
39# creates copy of code in output dir
40copy_code_files(project_slide.outputtimedir,__file__,dirname(project_slide.__file__)+sep+ project_slide.__name__+'.py' )
41myid = 0
42numprocs = 1
43start_screen_catcher(project_slide.outputtimedir, myid, numprocs)
44
45print 'USER:    ', project_slide.user
46#-------------------------------------------------------------------------------
47# Preparation of topographic data
48#
49# Convert ASC 2 DEM 2 PTS using source data and store result in source data
50#-------------------------------------------------------------------------------
51"""
52# filenames
53on_offshore10_dem_name = project_slide.on_offshore10_dem_name
54nsw_dem_name = project_slide.nsw_dem_name
55
56
57# creates DEM from asc data
58convert_dem_from_ascii2netcdf(on_offshore10_dem_name, use_cache=True, verbose=True)
59convert_dem_from_ascii2netcdf(nsw_dem_name, use_cache=True, verbose=True)
60
61#creates pts file for onshore DEM
62dem2pts(on_offshore10_dem_name, use_cache=True, verbose=True)
63dem2pts(nsw_dem_name,
64        easting_min=project_slide.eastingmin_nsw,
65        easting_max=project_slide.eastingmax_nsw,
66        northing_min=project_slide.northingmin_nsw,
67        northing_max= project_slide.northingmax_nsw,
68        use_cache=True, verbose=True)
69
70print 'create offshore'
71G11 = Geospatial_data(file_name = project_slide.offshore_dem_name1 + '.txt')
72G12 = Geospatial_data(file_name = project_slide.offshore_dem_name4 + '.txt')+\
73      Geospatial_data(file_name = project_slide.offshore_dem_name5 + '.txt')+\
74      Geospatial_data(file_name = project_slide.offshore_dem_name6 + '.txt')+\
75      Geospatial_data(file_name = project_slide.offshore_dem_name7 + '.txt')+\
76      Geospatial_data(file_name = project_slide.offshore_dem_name8 + '.txt')+\
77      Geospatial_data(file_name = project_slide.offshore_dem_name9 + '.txt')
78print 'create onshore'
79G2 = Geospatial_data(file_name = project_slide.on_offshore10_dem_name + '.pts')
80G4 = Geospatial_data(file_name = project_slide.nsw_dem_name + '.pts')
81print 'add'
82G = G11.clip(Geospatial_data(project_slide.poly_surveyclip)) +\
83    G12.clip(Geospatial_data(project_slide.polyAll)) +\
84    G2.clip(Geospatial_data(project_slide.poly_10mclip)) +\
85    (G4.clip(Geospatial_data(project_slide.polyAll))).clip_outside(Geospatial_data(project_slide.poly_surveyclip)).clip_outside(Geospatial_data(project_slide.poly_10mclip))
86print 'export points'
87G.export_points_file(project_slide.combined_dem_name + '.pts')
88#G.export_points_file(project_slide.combined_dem_name + '.xya')
89"""
90#----------------------------------------------------------------------------
91# Create the triangular mesh based on overall clipping polygon with a tagged
92# boundary and interior regions defined in project.py along with
93# resolutions (maximal area of per triangle) for each polygon
94#-------------------------------------------------------------------------------
95
96from anuga.pmesh.mesh_interface import create_mesh_from_regions
97meshname = project_slide.meshname+'.msh'
98remainder_res = 100000
99local_res = 25000
100gong_res = 500
101
102interior_regions = [[project_slide.poly_local, local_res],
103                    [project_slide.poly_gong, gong_res],
104                    [project_slide.poly_southgong, gong_res]]
105
106from caching import cache
107_ = cache(create_mesh_from_regions,
108          project_slide.polyAll,
109           {'boundary_tags': {'e0': [0], 'e1': [1], 'e2': [2],
110                              'e3': [3], 'e4': [4]},
111           'maximum_triangle_area': remainder_res,
112           'filename': meshname},
113           #'interior_regions': interior_regions},
114          verbose = True, evaluate=False)
115print 'created mesh'
116
117#-------------------------------------------------------------------------------                                 
118# Setup computational domain
119#-------------------------------------------------------------------------------                                 
120domain = Domain(meshname, use_cache = True, verbose = True)
121
122print 'Number of triangles = ', len(domain)
123print 'The extent is ', domain.get_extent()
124print domain.statistics()
125
126domain.set_name(project_slide.basename)
127domain.set_datadir(project_slide.outputtimedir)
128domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
129domain.set_minimum_storable_height(0.01)
130
131#-------------------------------------------------------------------------------                                 
132# Set elevation to mesh
133#-------------------------------------------------------------------------------
134
135tide = 0.0
136domain.set_quantity('elevation', 
137                    filename = project_slide.combined_dem_name + '.pts',
138                    use_cache = True,
139                    verbose = True,
140                    alpha = 0.1
141                    )
142
143#-------------------------------------------------------------------------------                                 
144# Setup boundary conditions
145#-------------------------------------------------------------------------------
146print 'Available boundary tags', domain.get_boundary_tags()
147
148Br = Reflective_boundary(domain)
149Bd = Dirichlet_boundary([tide,0,0])
150
151domain.set_boundary( {'e0': Bd,  'e1': Bd, 'e2': Bd, 'e3': Bd, 'e4': Bd} )
152
153#-------------------------------------------------------------------------------
154# Set up scenario (tsunami_source is a callable object used with set_quantity)
155#-------------------------------------------------------------------------------
156from smf import slide_tsunami
157
158# effect on a3D and wavelength
159gamma = 1.85
160massco = 1.0
161dragco = 1.0
162
163# no effect on a3D and wavelength but used in Double Gaussian
164#dx = 0.01
165kappa = 3.
166kappad = 0.8
167
168# this doesn't seem to apper anywhere in smf
169frictionco = 0.01
170
171# Bridgette to vary this parameter from 1 to 25 in steps of (?) by changing it in the project file
172                 
173tsunami_source = slide_tsunami(project_slide.bulli_length,
174                               project_slide.bulli_depth,
175                               project_slide.bulli_slope,
176                               width=project_slide.bulli_width,
177                               thickness=project_slide.bulli_thickness, 
178                               x0=project_slide.slide_origin_bulli[0], 
179                               y0=project_slide.slide_origin_bulli[1], 
180                               alpha=project_slide.bulli_alpha,
181                               dx=None,                     
182                               scale=project_slide.scale,
183                               domain=domain,
184                               verbose=True)
185
186
187#-------------------------------------------------------------------------------                                 
188# Setup initial conditions
189#-------------------------------------------------------------------------------
190
191domain.set_quantity('stage', tsunami_source)
192domain.set_quantity('friction', 0.0) 
193
194stage = domain.get_quantity('stage')
195etamax = stage.get_maximum_value()
196etamin = stage.get_minimum_value()
197
198filename = project_slide.outputtimedir + 'smf_data.csv'
199fid = open(filename,'w')
200s = 'scale, a3D, wavelength, etamax, etamin \n'
201fid.write(s)
202if project_slide.scale is not None:
203    s = '%.6f, %.6f, %.6f, %.6f, %.6f\n' \
204         %(project_slide.scale,tsunami_source.a3D,tsunami_source.wavelength,etamax,etamin)
205else:
206    s = '%s, %.6f, %.6f, %.6f, %.6f\n' \
207        %(project_slide.scale,tsunami_source.a3D,tsunami_source.wavelength,etamax,etamin)
208fid.write(s)   
209fid.close()
210
211#-------------------------------------------------------------------------------                                 
212# Evolve system through time
213#-------------------------------------------------------------------------------
214import time
215t0 = time.time()
216from Numeric import allclose
217from anuga.abstract_2d_finite_volumes.quantity import Quantity
218
219for t in domain.evolve(yieldstep = 30, finaltime = 5000): 
220    domain.write_time()
221       
222       
223print 'That took %.2f seconds' %(time.time()-t0)
224
225print 'finished'
Note: See TracBrowser for help on using the repository browser.