source: anuga_work/production/west_tas_2008/run_west_tas_smf.py @ 5442

Last change on this file since 5442 was 5442, checked in by ole, 15 years ago

Retired h-limiter and beta_h as per ticket:194.
All unit tests and validation tests pass.

File size: 7.1 KB
Line 
1"""Script for running a tsunami inundation scenario for Sydney, NSW, Australia.
2
3Source data such as elevation and boundary data is assumed to be available in
4directories specified by project_smf.py
5The output sww file is stored in project_smf.outputtimedir
6
7The scenario is defined by a triangular mesh created from project_smf.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
30from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
31
32# Application specific imports
33import project_smf              # Definition of file names and polygons
34
35#-------------------------------------------------------------------------------
36# Copy scripts to time stamped output directory and capture screen
37# output to file
38#-------------------------------------------------------------------------------
39
40# creates copy of code in output dir
41#copy_code_files(project_smf.outputtimedir,__file__,dirname(project_smf.__file__)+sep+ project_smf.__name__+'.py' )
42myid = 0
43numprocs = 1
44#start_screen_catcher(project_smf.outputtimedir, myid, numprocs)
45#barrier()
46
47print 'USER:    ', project_smf.user
48
49#-------------------------------------------------------------------------------
50# Preparation of topographic data
51#
52# Convert ASC 2 DEM 2 PTS using source data and store result in source data
53#-------------------------------------------------------------------------------
54
55# filenames
56onshore_250_dem_name = project_smf.onshore_250_dem_name
57meshname = project_smf.meshname+'.msh'
58
59# creates DEM from asc data
60convert_dem_from_ascii2netcdf(onshore_250_dem_name, use_cache=True, verbose=True)
61
62#creates pts file for onshore DEM
63dem2pts(onshore_250_dem_name, use_cache=True, verbose=True)
64
65print 'create offshore'
66G1 = Geospatial_data(file_name = project_smf.offshore_dem_name1 + '.txt')+\
67     Geospatial_data(file_name = project_smf.offshore_dem_name2 + '.txt')+\
68     Geospatial_data(file_name = project_smf.offshore_dem_name3 + '.txt')
69
70print 'create onshore'
71G2 = Geospatial_data(file_name = project_smf.onshore_250_dem_name + '.pts')
72
73print 'create coastline'
74G3 = Geospatial_data(file_name = project_smf.coast_line + '.txt')
75
76print 'add'
77G = G1.clip(Geospatial_data(project_smf.polyAll)) + G2 + G3
78
79print 'export points'
80G.export_points_file(project_smf.combined_dem_name + '.pts')
81
82
83#----------------------------------------------------------------------------
84# Create the triangular mesh based on overall clipping polygon with a tagged
85# boundary and interior regions defined in project_smf.py along with
86# resolutions (maximal area of per triangle) for each polygon
87#-------------------------------------------------------------------------------
88
89from anuga.pmesh.mesh_interface import create_mesh_from_regions
90remainder_res = 250000.
91region_res = 50000.
92local_res = 500.
93interior_regions = [[project_smf.poly_region, region_res],
94                    [project_smf.poly_local, local_res]]
95
96from caching import cache
97_ = cache(create_mesh_from_regions,
98          project_smf.polyAll,
99           {'boundary_tags': {'e0': [0], 'e1': [1], 'e2': [2],
100                              'e3': [3], 'e4':[4]},
101           'maximum_triangle_area': remainder_res,
102           'filename': meshname,
103           'interior_regions': interior_regions},
104          verbose = True, evaluate=False)
105print 'created mesh'
106
107#-------------------------------------------------------------------------------                                 
108# Setup computational domain
109#-------------------------------------------------------------------------------                                 
110domain = Domain(meshname, use_cache = True, verbose = True)
111
112print 'Number of triangles = ', len(domain)
113print 'The extent is ', domain.get_extent()
114print domain.statistics()
115
116domain.set_name(project_smf.basename)
117domain.set_datadir(project_smf.outputtimedir)
118domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
119domain.set_minimum_storable_height(0.01)
120domain.set_store_vertices_uniquely(False)
121
122
123#-------------------------------------------------------------------------------                                 
124# Setup initial conditions
125#-------------------------------------------------------------------------------
126
127tide = 0.0
128domain.set_quantity('stage', tide)
129domain.set_quantity('friction', 0.0) 
130domain.set_quantity('elevation', 
131                    filename = project_smf.combined_dem_name + '.pts',
132                    use_cache = True,
133                    verbose = True,
134                    alpha = 0.1
135                    )
136
137#-------------------------------------------------------------------------------
138# Set up scenario (tsunami_source is a callable object used with set_quantity)
139#-------------------------------------------------------------------------------
140from smf import slide_tsunami
141
142tsunami_source = slide_tsunami(length=project_smf.length,
143                               width=project_smf.width,
144                               depth=project_smf.depth,
145                               slope=project_smf.slope,
146                               thickness=project_smf.thickness, 
147                               x0=project_smf.smf_origin[0], 
148                               y0=project_smf.smf_origin[1], 
149                               alpha=project_smf.alpha, 
150                               domain=domain)
151
152#-------------------------------------------------------------------------------                                 
153# Setup boundary conditions
154#-------------------------------------------------------------------------------
155print 'Available boundary tags', domain.get_boundary_tags()
156
157Br = Reflective_boundary(domain)
158Bd = Dirichlet_boundary([tide,0,0])
159
160domain.set_boundary( {'e0': Bd,  'e1': Bd, 'e2': Bd, 'e3': Bd, 'e4': Bd} )
161
162
163#-------------------------------------------------------------------------------                                 
164# Evolve system through time
165#-------------------------------------------------------------------------------
166import time
167from Numeric import allclose
168from anuga.abstract_2d_finite_volumes.quantity import Quantity
169
170t0 = time.time()
171
172for t in domain.evolve(yieldstep = 30, finaltime = 5000): 
173    domain.write_time()
174    domain.write_boundary_statistics(tags = 'e2')
175    stagestep = domain.get_quantity('stage') 
176
177    if allclose(t, 30):
178        smf = Quantity(domain)
179        smf.set_values(tsunami_source)
180        domain.set_quantity('stage', smf + stagestep)
181   
182print 'That took %.2f seconds' %(time.time()-t0)
183
184print 'finished'
Note: See TracBrowser for help on using the repository browser.