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

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

Bridgette's study for investigating sensitivity of input wave at Wollongong

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# Setup initial conditions
188#-------------------------------------------------------------------------------
189
190domain.set_quantity('stage', tsunami_source)
191domain.set_quantity('friction', 0.0) 
192
193stage = domain.get_quantity('stage')
194etamax = stage.get_maximum_value()
195etamin = stage.get_minimum_value()
196
197filename = project_slide.outputtimedir + 'smf_data.csv'
198fid = open(filename,'w')
199s = 'scale, a3D, wavelength, etamax, etamin \n'
200fid.write(s)
201if project_slide.scale is not None:
202    s = '%.6f, %.6f, %.6f, %.6f, %.6f\n' \
203         %(project_slide.scale,tsunami_source.a3D,tsunami_source.wavelength,etamax,etamin)
204else:
205    s = '%s, %.6f, %.6f, %.6f, %.6f\n' \
206        %(project_slide.scale,tsunami_source.a3D,tsunami_source.wavelength,etamax,etamin)
207fid.write(s)   
208fid.close()
209
210#-------------------------------------------------------------------------------                                 
211# Evolve system through time
212#-------------------------------------------------------------------------------
213import time
214t0 = time.time()
215from Numeric import allclose
216from anuga.abstract_2d_finite_volumes.quantity import Quantity
217
218for t in domain.evolve(yieldstep = 30, finaltime = 5000): 
219    domain.write_time()
220       
221       
222print 'That took %.2f seconds' %(time.time()-t0)
223
224print 'finished'
Note: See TracBrowser for help on using the repository browser.