source: anuga_work/production/gold_coast_2007/run_goldcoast.py @ 4347

Last change on this file since 4347 was 4310, checked in by sexton, 18 years ago

update damage scripts for FESA areas

File size: 5.9 KB
Line 
1"""Script for running a hypothetical inundation scenario for Gold Coast,
2QLD, Australia.
3
4Source data such as elevation and boundary data is assumed to be available in
5directories specified by project.py
6The output sww file is stored in project.outputtimedir
7
8The scenario is defined by a triangular mesh created from project.polygon,
9the elevation data and a tsunami wave generated by s submarine mass failure.
10
11Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis, GA - 2006
12"""
13
14#-------------------------------------------------------------------------------
15# Import necessary modules
16#-------------------------------------------------------------------------------
17
18# Standard modules
19import os
20import time
21from shutil import copy
22from os.path import dirname, basename
23from os import mkdir, access, F_OK, sep
24import sys
25
26# Related major packages
27from anuga.shallow_water import Domain, Reflective_boundary, Dirichlet_boundary, Time_boundary
28from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
29from anuga.geospatial_data.geospatial_data import *
30from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
31
32# Application specific imports
33import project              # 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
41copy_code_files(project.outputtimedir,__file__,dirname(project.__file__)+sep+ project.__name__+'.py' )
42myid = 0
43numprocs = 1
44start_screen_catcher(project.outputtimedir, myid, numprocs)
45
46print 'USER:    ', project.user
47
48#-------------------------------------------------------------------------------
49# Preparation of topographic data
50#
51# Convert ASC 2 DEM 2 PTS using source data and store result in source data
52#-------------------------------------------------------------------------------
53
54# filenames
55gc_dem_name = project.gc_dem_name
56meshname = project.meshname+'.msh'
57'''
58# creates DEM from asc data
59convert_dem_from_ascii2netcdf(gc_dem_name, use_cache=True, verbose=True)
60
61#creates pts file for onshore DEM
62dem2pts(gc_dem_name, use_cache=True, verbose=True)
63
64G = Geospatial_data(file_name = project.gc_dem_name + '.pts') + \
65    Geospatial_data(file_name = project.offshore_name + '.xya') +\
66    Geospatial_data(file_name = project.coast_name + '.xya')
67
68G.export_points_file(project.combined_dem_name + '.pts')
69'''
70#----------------------------------------------------------------------------
71# Create the triangular mesh based on overall clipping polygon with a tagged
72# boundary and interior regions defined in project.py along with
73# resolutions (maximal area of per triangle) for each polygon
74#-------------------------------------------------------------------------------
75
76from anuga.pmesh.mesh_interface import create_mesh_from_regions
77remainder_res = 1000000.
78int_res = 10000.
79interior_regions = [[project.poly_int, int_res]]
80
81from caching import cache
82_ = cache(create_mesh_from_regions,
83          project.polyAll,
84           {'boundary_tags': {'e0': [0], 'e1': [1], 'e2': [2], 'e3': [3]},
85           'maximum_triangle_area': remainder_res,
86           'filename': meshname,
87           'interior_regions': interior_regions},
88          verbose = True, evaluate=False)
89
90#-------------------------------------------------------------------------------                                 
91# Setup computational domain
92#-------------------------------------------------------------------------------                                 
93domain = Domain(meshname, use_cache = True, verbose = True)
94
95print 'Number of triangles = ', len(domain)
96print 'The extent is ', domain.get_extent()
97print domain.statistics()
98
99domain.set_name(project.basename)
100domain.set_datadir(project.outputtimedir)
101domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
102domain.set_minimum_storable_height(0.01)
103
104#-------------------------------------------------------------------------------                                 
105# Setup initial conditions
106#-------------------------------------------------------------------------------
107
108#tide = 3.0
109tide = 1.0
110from polygon import *
111#IC = Polygon_function( [(project.poly_bathy, 0.)], default = tide)
112#IC = Polygon_function( [(poly, 0.)], default = tide)
113from anuga.utilities.polygon import read_polygon
114poly = read_polygon('test.csv')
115IC = Polygon_function( [(poly, 1000.)], default = 0.0)
116domain.set_quantity('stage', IC)
117domain.set_quantity('friction', 0.0) 
118domain.set_quantity('elevation', 
119                    filename = project.combined_dem_name + '.pts',
120                    use_cache = True,
121                    verbose = True,
122                    alpha = 0.1
123                    )
124
125#-------------------------------------------------------------------------------                                 
126# Setup boundary conditions
127#-------------------------------------------------------------------------------
128print 'Available boundary tags', domain.get_boundary_tags()
129
130Br = Reflective_boundary(domain)
131Bd = Dirichlet_boundary([tide,0,0])
132# 10 min period = 1/600 frequency
133from math import sin
134Bw = Time_boundary(domain=domain, # Time dependent boundary
135                   #f=lambda t: [1.0*sin(t/600.), 0, 0])
136                   f=lambda t: [(60<t<3660)*0.1, 0, 0])
137
138domain.set_boundary( {'e0': Bw,  'e1': Bd, 'e2': Bd, 'e3': Bd} )
139
140
141#-------------------------------------------------------------------------------                                 
142# Evolve system through time
143#-------------------------------------------------------------------------------
144import time
145
146t0 = time.time()
147
148for t in domain.evolve(yieldstep = 10, finaltime = 15000): 
149    domain.write_time()
150    domain.write_boundary_statistics(tags = 'e0')
151   
152print 'That took %.2f seconds' %(time.time()-t0)
153
154print 'finished'
Note: See TracBrowser for help on using the repository browser.