source: anuga_work/production/wollongong_2006/run_kembla.py @ 4696

Last change on this file since 4696 was 4696, checked in by sexton, 17 years ago

tests for Lex

File size: 5.1 KB
Line 
1"""Script for running a tsunami inundation scenario for Port Kembla, NSW, Australia.
2
3Source data such as elevation and boundary data is assumed to be available in
4directories specified by project_kembla.py
5
6The scenario is defined by a triangular mesh created from project_kembla.polygon,
7the elevation data and an input wave.
8
9Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis and Jane Sexton, GA - 2006
10"""
11
12#-------------------------------------------------------------------------------
13# Import necessary modules
14#-------------------------------------------------------------------------------
15
16# Standard modules
17import os
18import time
19from shutil import copy
20from os.path import dirname, basename
21from os import mkdir, access, F_OK, sep
22import sys
23
24# Related major packages
25from anuga.shallow_water import Domain, Reflective_boundary, Dirichlet_boundary, Time_boundary
26from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
27from anuga.geospatial_data.geospatial_data import *
28from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
29
30# Application specific imports
31import project_kembla              # Definition of file names and polygons
32
33#-------------------------------------------------------------------------------
34# Preparation of topographic data
35#
36# Convert ASC 2 DEM 2 PTS using source data and store result in source data
37#-------------------------------------------------------------------------------
38
39# filenames
40dem_name = project_kembla.dem_name
41meshname = project_kembla.meshname+'.msh'
42
43# creates DEM from asc data
44convert_dem_from_ascii2netcdf(dem_name, use_cache=True, verbose=True)
45
46#creates pts file for onshore DEM
47dem2pts(dem_name, use_cache=True, verbose=True)
48
49print 'create onshore'
50G = Geospatial_data(file_name = project_kembla.dem_name + '.pts')
51
52print 'export points'
53G.export_points_file(project_kembla.combined_dem_name + '.pts')
54#G.export_points_file(project_kembla.combined_dem_name + '.xya')
55
56#----------------------------------------------------------------------------
57# Create the triangular mesh based on overall clipping polygon with a tagged
58# boundary and interior regions defined in project.py along with
59# resolutions (maximal area of per triangle) for each polygon
60#-------------------------------------------------------------------------------
61
62from anuga.pmesh.mesh_interface import create_mesh_from_regions
63remainder_res = 100000
64bay_res = 5000
65interior_regions = [[project_kembla.poly_bay, bay_res]]
66
67from caching import cache
68_ = cache(create_mesh_from_regions,
69          project_kembla.polyAll,
70           {'boundary_tags': {'top': [0], 'right': [1], 'bottom': [2],
71                              'left': [3]},
72           'maximum_triangle_area': remainder_res,
73           'filename': meshname,
74           'interior_regions': interior_regions},
75          verbose = True, evaluate=False)
76print 'created mesh'
77
78#-------------------------------------------------------------------------------                                 
79# Setup computational domain
80#-------------------------------------------------------------------------------                                 
81domain = Domain(meshname, use_cache = True, verbose = True)
82
83print 'Number of triangles = ', len(domain)
84print 'The extent is ', domain.get_extent()
85print domain.statistics()
86
87domain.set_name(project_kembla.basename)
88domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
89domain.set_minimum_storable_height(0.01)
90
91#-------------------------------------------------------------------------------                                 
92# Setup initial conditions
93#-------------------------------------------------------------------------------
94
95tide = 0.0
96domain.set_quantity('stage', tide)
97domain.set_quantity('friction', 0.0) 
98domain.set_quantity('elevation', 
99                    filename = project_kembla.combined_dem_name + '.pts',
100                    use_cache = True,
101                    verbose = True,
102                    alpha = 0.1
103                    )
104
105#-------------------------------------------------------------------------------                                 
106# Setup boundary conditions
107#-------------------------------------------------------------------------------
108print 'Available boundary tags', domain.get_boundary_tags()
109
110Br = Reflective_boundary(domain)
111Bd = Dirichlet_boundary([tide,0,0])
112# 10 min square wave starting at 1 min, 6m high
113Bw = Time_boundary(domain=domain,
114                   f=lambda t: [(60<t<660)*6, 0, 0])
115
116domain.set_boundary( {'top': Bd,  'right': Bw, 'bottom': Bd, 'left': Bd} )
117
118
119#-------------------------------------------------------------------------------                                 
120# Evolve system through time
121#-------------------------------------------------------------------------------
122import time
123t0 = time.time()
124from Numeric import allclose
125from anuga.abstract_2d_finite_volumes.quantity import Quantity
126
127for t in domain.evolve(yieldstep = 30, finaltime = 5000): 
128    domain.write_time()
129    domain.write_boundary_statistics(tags = '')
130    stagestep = domain.get_quantity('stage') 
131   
132print 'That took %.2f seconds' %(time.time()-t0)
133
134print 'finished'
Note: See TracBrowser for help on using the repository browser.