source: production/karratha_2005/run_karratha.py @ 1855

Last change on this file since 1855 was 1855, checked in by ole, 19 years ago

Better interface to meshing for karratha + some pypar work

File size: 4.6 KB
RevLine 
[1845]1"""Script for running a tsunami inundation scenario for Karratha, WA, Australia.
2
3Source data such as elevation and boundary data is assumed to be available in
4directories specified by project.py
5The output sww file is stored in project.outputdir
6
7The scenario is defined by a triangular mesh created from project.polygon,
8the elevation data and boundary data obtained from a tsunami simulation done with MOST.
9
10Ole Nielsen, GA - 2005
[1781]11"""
12
[1855]13tide = 0.75   #HMWS estimate by Colin French, GA
14
15
[1781]16import os
17import time
18
[1784]19
[1845]20from pyvolution.shallow_water import Domain, Reflective_boundary, File_boundary,\
21     Dirichlet_boundary, Time_boundary, Transmissive_boundary
[1796]22from pyvolution.data_manager import convert_dem_from_ascii2netcdf,\
23     dem2pts, ferret2sww
[1784]24from pyvolution.pmesh2domain import pmesh_to_domain_instance
25from caching import cache
[1781]26import project
27
[1796]28#Data preparation
[1781]29#Convert ASC 2 DEM 2 PTS using source data and store result in source data
[1784]30demname = project.demname
31
32cache(convert_dem_from_ascii2netcdf, demname, {'verbose': True},
33      dependencies = [demname + '.asc'],
34      verbose = True)
[1796]35      #evaluate = True)
[1784]36
37cache(dem2pts, demname, {'verbose': True},
38      dependencies = [demname + '.dem'],     
39      verbose = True)
40
[1801]41
[1796]42#Convert MOST boundary
[1801]43source_dir = project.boundarydir
[1784]44
[1801]45from pyvolution.data_manager import ferret2sww
[1796]46
[1801]47south = project.south
48north = project.north
49west = project.west
50east = project.east
51
52
[1855]53cache(ferret2sww,
54      (source_dir+project.boundary_basename,
55       project.boundary_basename), 
56      {'verbose': True,
57       'minlat': south-1,
58       'maxlat': north+1,
59       'minlon': west-1,
60       'maxlon': east+1,
61       'origin': project.mesh_origin,
62       'mean_stage': tide,
63       'zscale': 1,
64       'fail_on_NaN': False,
65       'inverted_bathymetry': True},
66      verbose = True)
67#FIXME: Dependencies
[1845]68
[1855]69   
[1843]70#ferret2sww(source_dir+project.boundary_basename,
71#           project.boundary_basename,
72#           verbose=True,
73#           minlat=south-1, maxlat=north+1,
74#           minlon=west-1, maxlon=east+1,
75#           origin = project.mesh_origin,
[1855]76#           mean_stage = tide,
[1843]77#           zscale = 1,
78#           fail_on_NaN = False,
[1855]79#           inverted_bathymetry = True)
[1832]80
[1801]81
[1855]82#Create Triangular Mesh
83from create_mesh import create_mesh
[1801]84
[1855]85interior_regions = [[project.karratha_polygon, 10000]]
86m = cache(create_mesh,
87          project.polygon,
88          {'boundary_tags': {'back': [7, 8], 'side': [0, 6], 'ocean': [1, 2, 3, 4, 5]},
89          'resolution': 80000,
90          'filename': project.meshname + '.msh',
91          'interior_regions': interior_regions},
92          verbose = True)         
93          #verbose = True,
94          #evaluate = True )
95
96
[1801]97#Setup domain
[1784]98mesh = project.meshname + '.msh'
99domain = cache(pmesh_to_domain_instance, (mesh, Domain),
100               dependencies = [mesh],                     
101               verbose = True)               
102
103
104domain.set_name(project.basename)
[1786]105domain.set_datadir(project.outputdir)
106domain.store = True
107
[1855]108#domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
109domain.quantities_to_be_stored = ['stage']
110       
[1832]111print 'Number of triangles = ', len(domain)
[1785]112print 'The extent is ', domain.get_extent()
[1784]113
[1796]114
[1845]115
116#Setup Initial Conditions
[1843]117domain.set_quantity('friction', 0)
[1832]118domain.set_quantity('stage', tide)
[1786]119domain.set_quantity('elevation',
120                    filename = demname + '.pts',
121                    use_cache = True,
122                    verbose = True)
[1784]123
[1801]124
[1845]125
126#Setup Boundary Conditions
[1830]127print domain.get_boundary_tags()
128
[1832]129Bf = File_boundary(project.boundary_basename + '.sww', domain, verbose = True)
[1855]130#domain.starttime = 3000  #Obtained from MOST
131domain.starttime = 0  #Obtained from MOST
[1830]132
[1784]133Br = Reflective_boundary(domain)
[1832]134Bt = Transmissive_boundary(domain)
135Bd = Dirichlet_boundary([tide,0,0])
[1830]136Bw = Time_boundary(domain=domain,
[1832]137                   f=lambda t: [(60<t<660)*4, 0, 0])
[1784]138
[1832]139domain.set_boundary( {'back': Br,'side': Bd, 'ocean': Bf} )
[1784]140
[1843]141#domain.set_boundary( {'back': Br,'side': Bd, 'ocean': Bd} )  #Nothing
[1832]142
[1843]143
[1784]144#Run
[1855]145#for t in domain.evolve(yieldstep = 600, finaltime = 15000):
146#    domain.write_time()
147#    domain.write_boundary_statistics(tags = 'ocean') #quantities = 'stage')
148#
149#for t in domain.evolve(yieldstep = 10, finaltime = 35000):
150#    domain.write_time()
151#    domain.write_boundary_statistics(tags = 'ocean') #quantities = 'stage')
152
[1843]153for t in domain.evolve(yieldstep = 60, finaltime = 40000):
[1784]154    domain.write_time()
[1855]155    domain.write_boundary_statistics(tags = 'ocean') #quantities = 'stage')       
Note: See TracBrowser for help on using the repository browser.