source: trunk/anuga_work/development/gong_2008/run_gong_slide.py @ 7884

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

updates for gong study

File size: 9.5 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
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
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
40#-------------------------------------------------------------------------------
41# Preparation of topographic data
42#
43# Convert ASC 2 DEM 2 PTS using source data and store result in source data
44#-------------------------------------------------------------------------------
45"""
46# filenames
47on_offshore10_dem_name = project_slide.on_offshore10_dem_name
48nsw_dem_name = project_slide.nsw_dem_name
49
50
51# creates DEM from asc data
52convert_dem_from_ascii2netcdf(on_offshore10_dem_name, use_cache=True, verbose=True)
53convert_dem_from_ascii2netcdf(nsw_dem_name, use_cache=True, verbose=True)
54
55#creates pts file for onshore DEM
56dem2pts(on_offshore10_dem_name, use_cache=True, verbose=True)
57dem2pts(nsw_dem_name,
58        easting_min=project_slide.eastingmin_nsw,
59        easting_max=project_slide.eastingmax_nsw,
60        northing_min=project_slide.northingmin_nsw,
61        northing_max= project_slide.northingmax_nsw,
62        use_cache=True, verbose=True)
63
64print 'create offshore'
65G11 = Geospatial_data(file_name = project_slide.offshore_dem_name1 + '.txt')
66G12 = Geospatial_data(file_name = project_slide.offshore_dem_name4 + '.txt')+\
67      Geospatial_data(file_name = project_slide.offshore_dem_name5 + '.txt')+\
68      Geospatial_data(file_name = project_slide.offshore_dem_name6 + '.txt')+\
69      Geospatial_data(file_name = project_slide.offshore_dem_name7 + '.txt')+\
70      Geospatial_data(file_name = project_slide.offshore_dem_name8 + '.txt')+\
71      Geospatial_data(file_name = project_slide.offshore_dem_name9 + '.txt')
72print 'create onshore'
73G2 = Geospatial_data(file_name = project_slide.on_offshore10_dem_name + '.pts')
74G4 = Geospatial_data(file_name = project_slide.nsw_dem_name + '.pts')
75print 'add'
76G = G11.clip(Geospatial_data(project_slide.poly_surveyclip)) +\
77    G12.clip(Geospatial_data(project_slide.polyAll)) +\
78    G2.clip(Geospatial_data(project_slide.poly_10mclip)) +\
79    (G4.clip(Geospatial_data(project_slide.polyAll))).clip_outside(Geospatial_data(project_slide.poly_surveyclip)).clip_outside(Geospatial_data(project_slide.poly_10mclip))
80print 'export points'
81G.export_points_file(project_slide.combined_dem_name + '.pts')
82#G.export_points_file(project_slide.combined_dem_name + '.xya')
83"""
84#----------------------------------------------------------------------------
85# Create the triangular mesh based on overall clipping polygon with a tagged
86# boundary and interior regions defined in project.py along with
87# resolutions (maximal area of per triangle) for each polygon
88#-------------------------------------------------------------------------------
89
90from anuga.pmesh.mesh_interface import create_mesh_from_regions
91meshname = project_slide.meshname+'.msh'
92remainder_res = 100000*100
93local_res = 25000*1000
94gong_res = 500*1000
95interior_regions = [[project_slide.poly_local, local_res],
96                    [project_slide.poly_gong, gong_res],
97                    [project_slide.poly_southgong, gong_res]]
98
99from caching import cache
100_ = cache(create_mesh_from_regions,
101          project_slide.polyAll,
102           {'boundary_tags': {'e0': [0], 'e1': [1], 'e2': [2],
103                              'e3': [3], 'e4':[4]},
104           'maximum_triangle_area': remainder_res,
105           'filename': meshname,
106           'interior_regions': interior_regions},
107          verbose = True, evaluate=False)
108print 'created mesh'
109
110#-------------------------------------------------------------------------------                                 
111# Setup computational domain
112#-------------------------------------------------------------------------------                                 
113domain = Domain(meshname, use_cache = True, verbose = True)
114
115print 'Number of triangles = ', len(domain)
116print 'The extent is ', domain.get_extent()
117print domain.statistics()
118
119domain.set_name(project_slide.basename)
120domain.set_datadir(project_slide.outputtimedir)
121domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
122domain.set_minimum_storable_height(0.01)
123
124#-------------------------------------------------------------------------------                                 
125# Set elevation to mesh
126#-------------------------------------------------------------------------------
127
128tide = 0.0
129domain.set_quantity('elevation', 
130                    filename = project_slide.combined_dem_name + '.pts',
131                    use_cache = True,
132                    verbose = True,
133                    alpha = 0.1
134                    )
135
136#-------------------------------------------------------------------------------                                 
137# Setup boundary conditions
138#-------------------------------------------------------------------------------
139print 'Available boundary tags', domain.get_boundary_tags()
140
141Br = Reflective_boundary(domain)
142Bd = Dirichlet_boundary([tide,0,0])
143
144domain.set_boundary( {'e0': Bd,  'e1': Bd, 'e2': Bd, 'e3': Bd, 'e4': Bd} )
145
146
147#-------------------------------------------------------------------------------
148# Set up scenario (tsunami_source is a callable object used with set_quantity)
149#-------------------------------------------------------------------------------
150from smf import slide_tsunami
151
152# effect on a3D and wavelength
153l0 = project_slide.birubi_length
154w0 = project_slide.birubi_width
155d0 = project_slide.birubi_depth
156s0 = project_slide.birubi_slope
157t0 = project_slide.birubi_thickness
158a0 = project_slide.birubi_alpha
159x0 = project_slide.slide_origin_birubi_a[0]
160y0 = project_slide.slide_origin_birubi_a[1]
161length = l0
162width = w0
163depth = d0
164slope = s0
165thickness = t0
166alpha = a0
167gamma0 = 1.85
168gamma = gamma0
169m0 = 1.
170massco = m0
171drag0 = 1.
172dragco = drag0
173
174# no effect on a3D and wavelength but used in Double Gaussian
175dx = 0.01
176kappa0 = 3.
177kappad0 = 0.8
178kappa = kappa0
179kappad = kappad0
180
181# this doesn't seem to apper anywhere in smf
182frictionc0 = 0.01
183frictionco = frictionc0
184
185# scaling for Double Gaussian function
186scale0 = 100. # Bridgette's fiddle
187scale = scale0 + 10.
188
189for i in range(1):
190    scale = scale - 10.   
191    mydir = project_slide.outputdir+'testing_' + str(int(scale))
192    print 'dir_comment', mydir
193    start_screen_catcher(mydir, 0, 1)
194    # creates copy of code in output dir
195    copy_code_files(mydir,__file__,dirname(project_slide.__file__)+sep+ project_slide.__name__+'.py' )
196    start_screen_catcher(mydir, 0, 1)
197
198    print 'USER:    ', project_slide.user
199
200    tsunami_source = slide_tsunami(length=length,
201                                   width=width,
202                                   depth=depth,
203                                   slope=slope,
204                                   thickness=thickness, 
205                                   x0=x0, 
206                                   y0=y0, 
207                                   alpha=alpha,
208                                   gamma=gamma,
209                                   massco=massco,
210                                   dragco=dragco,
211                                   frictionco=frictionco,
212                                   dx=dx,
213                                   kappa=kappa,
214                                   kappad=kappad,
215                                   scale=scale,
216                                   domain=domain,
217                                   verbose=True)
218
219    print 'hello', scale, tsunami_source.wavelength, tsunami_source.a3D
220   
221    #-------------------------------------------------------------------------------                                 
222    # Setup initial conditions
223    #-------------------------------------------------------------------------------
224
225    domain.set_quantity('stage', tsunami_source)
226    domain.set_quantity('friction', 0.0) 
227
228
229    #-------------------------------------------------------------------------------                                 
230    # Evolve system through time
231    #-------------------------------------------------------------------------------
232    import time
233    t0 = time.time()
234    from Numeric import allclose
235    from anuga.abstract_2d_finite_volumes.quantity import Quantity
236
237    for t in domain.evolve(yieldstep = 30, finaltime = 5000): 
238        domain.write_time()
239        domain.write_boundary_statistics(tags = 'e14')
240       
241       
242    print 'That took %.2f seconds' %(time.time()-t0)
243
244    print 'finished'
245
Note: See TracBrowser for help on using the repository browser.