[3559] | 1 | """Script for running a tsunami inundation scenario for Hobart, TAS, Australia. |
---|
| 2 | |
---|
| 3 | Source data such as elevation and boundary data is assumed to be available in |
---|
| 4 | directories specified by project.py |
---|
| 5 | The output sww file is stored in project.outputtimedir |
---|
| 6 | |
---|
| 7 | The scenario is defined by a triangular mesh created from project.polygon, |
---|
[3615] | 8 | the elevation data and a tsunami wave generated by MOST. |
---|
[3559] | 9 | |
---|
| 10 | Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis, GA - 2006 |
---|
| 11 | """ |
---|
| 12 | #-------------------------------------------------------------------------------# Import necessary modules |
---|
| 13 | #------------------------------------------------------------------------------- |
---|
| 14 | |
---|
| 15 | # Standard modules |
---|
| 16 | import os |
---|
| 17 | import time |
---|
| 18 | from shutil import copy |
---|
| 19 | from os import mkdir, access, F_OK |
---|
| 20 | import sys |
---|
| 21 | |
---|
| 22 | # Related major packages |
---|
[3626] | 23 | from anuga.shallow_water import Domain, Reflective_boundary, \ |
---|
[3559] | 24 | Dirichlet_boundary, Time_boundary, File_boundary |
---|
[3626] | 25 | from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts |
---|
[3650] | 26 | from anuga.abstract_2d_finite_volumes.combine_pts import combine_rectangular_points_files |
---|
[3559] | 27 | from anuga.geospatial_data.geospatial_data import * |
---|
[3626] | 28 | from anuga.abstract_2d_finite_volumes.util import Screen_Catcher |
---|
[3559] | 29 | |
---|
| 30 | # Application specific imports |
---|
| 31 | import project # Definition of file names and polygons |
---|
| 32 | |
---|
| 33 | #------------------------------------------------------------------------------- |
---|
| 34 | # Copy scripts to time stamped output directory and capture screen |
---|
| 35 | # output to file |
---|
| 36 | #------------------------------------------------------------------------------- |
---|
| 37 | |
---|
| 38 | # creates copy of code in output dir if dir doesn't exist |
---|
| 39 | if access(project.outputtimedir,F_OK) == 0 : |
---|
| 40 | mkdir (project.outputtimedir) |
---|
| 41 | copy (project.codedirname, project.outputtimedir + project.codename) |
---|
[3626] | 42 | copy (project.codedir + 'run_hobart.py', project.outputtimedir + 'run_hobart.py') |
---|
[3559] | 43 | print'output dir', project.outputtimedir |
---|
| 44 | |
---|
| 45 | #normal screen output is stored in |
---|
| 46 | screen_output_name = project.outputtimedir + "screen_output.txt" |
---|
| 47 | screen_error_name = project.outputtimedir + "screen_error.txt" |
---|
| 48 | |
---|
| 49 | #used to catch screen output to file |
---|
| 50 | sys.stdout = Screen_Catcher(screen_output_name) |
---|
| 51 | #sys.stderr = Screen_Catcher(screen_output_name) |
---|
| 52 | sys.stderr = Screen_Catcher(screen_error_name) |
---|
| 53 | |
---|
| 54 | print 'USER: ', project.user |
---|
| 55 | |
---|
| 56 | #------------------------------------------------------------------------------- |
---|
| 57 | # Preparation of topographic data |
---|
| 58 | # |
---|
| 59 | # Convert ASC 2 DEM 2 PTS using source data and store result in source data |
---|
| 60 | #------------------------------------------------------------------------------- |
---|
| 61 | |
---|
| 62 | # filenames |
---|
[3671] | 63 | onshore_offshore_dem_name = project.onshore_offshore_dem_name |
---|
[3683] | 64 | onshore_offshore_dem_name_25 = project.onshore_offshore_dem_name_25 |
---|
[3559] | 65 | meshname = project.meshname+'.msh' |
---|
| 66 | source_dir = project.boundarydir |
---|
| 67 | |
---|
| 68 | copied_files = False |
---|
| 69 | |
---|
[3683] | 70 | # create DEM from 50m asc data |
---|
[3671] | 71 | convert_dem_from_ascii2netcdf(onshore_offshore_dem_name, use_cache=True, verbose=True) |
---|
[3660] | 72 | |
---|
[3683] | 73 | # creates pts file for combined 50m DEM |
---|
[3671] | 74 | dem2pts(onshore_offshore_dem_name, use_cache=True, verbose=True) |
---|
[3650] | 75 | |
---|
[3683] | 76 | # 25m data (clipping the around the Hobart area) |
---|
| 77 | convert_dem_from_ascii2netcdf(onshore_offshore_dem_name_25, use_cache=True, verbose=True) |
---|
| 78 | # creates pts file for 25m data around Hobart |
---|
| 79 | dem2pts(onshore_offshore_dem_name_25, project.hobart_dem_name_25, |
---|
| 80 | easting_min=project.eastingmin25, |
---|
| 81 | easting_max=project.eastingmax25, |
---|
| 82 | northing_min=project.northingmin25, |
---|
| 83 | northing_max= project.northingmax25, |
---|
| 84 | use_cache=True, |
---|
| 85 | verbose=True) |
---|
| 86 | |
---|
| 87 | # combining the 50m and Hobart 25m data |
---|
| 88 | combine_rectangular_points_files(project.hobart_dem_name_25 + '.pts', |
---|
| 89 | project.onshore_offshore_dem_name + '.pts', |
---|
| 90 | project.combined_dem_name + '.pts') |
---|
| 91 | |
---|
| 92 | # 25m data (clipping the around site 24 on Bruny Island) |
---|
| 93 | convert_dem_from_ascii2netcdf(onshore_offshore_dem_name_25, use_cache=True, verbose=True) |
---|
| 94 | # creates pts file for 25m data around site 24 at Bruny Island |
---|
| 95 | dem2pts(onshore_offshore_dem_name_25, project.bruny_dem_name_25, |
---|
| 96 | easting_min=project.eastingmin25_2, |
---|
| 97 | easting_max=project.eastingmax25_2, |
---|
| 98 | northing_min=project.northingmin25_2, |
---|
| 99 | northing_max= project.northingmax25_2, |
---|
| 100 | use_cache=True, |
---|
| 101 | verbose=True) |
---|
| 102 | |
---|
| 103 | # combining the 50m and Hobart 25m data with Bruny Island 25m data |
---|
| 104 | combine_rectangular_points_files(project.bruny_dem_name_25 + '.pts', |
---|
| 105 | project.combined_dem_name + '.pts', |
---|
| 106 | project.combined_dem_name_2 + '.pts') |
---|
| 107 | |
---|
[3671] | 108 | # create geospatial data set and export |
---|
[3683] | 109 | #G = Geospatial_data(file_name = project.onshore_offshore_dem_name + '.pts') |
---|
| 110 | #G.export_points_file(project.combined_dem_name + '.pts') |
---|
[3661] | 111 | |
---|
| 112 | #---------------------------------------------------------------------------- |
---|
[3559] | 113 | # Create the triangular mesh based on overall clipping polygon with a tagged |
---|
| 114 | # boundary and interior regions defined in project.py along with |
---|
| 115 | # resolutions (maximal area of per triangle) for each polygon |
---|
| 116 | #------------------------------------------------------------------------------- |
---|
| 117 | |
---|
| 118 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
| 119 | |
---|
[3615] | 120 | # use 75 for onshore components (12.5m DEM) |
---|
[3626] | 121 | island_res = 35000 |
---|
[3695] | 122 | hobart_res = 5000 |
---|
[3626] | 123 | peninsula_res = 35000 |
---|
[3679] | 124 | interior_regions = [[project.poly_hobart1, hobart_res], |
---|
| 125 | [project.poly_hobart2, hobart_res], |
---|
[3695] | 126 | [project.poly_hobart3, hobart_res], |
---|
| 127 | [project.poly_hobart4, hobart_res]] |
---|
[3559] | 128 | |
---|
| 129 | print 'number of interior regions', len(interior_regions) |
---|
| 130 | |
---|
| 131 | from caching import cache |
---|
| 132 | _ = cache(create_mesh_from_regions, |
---|
| 133 | project.polyAll, |
---|
[3673] | 134 | {'boundary_tags': {'e0': [0], 'e1': [1], 'e2': [2], |
---|
[3679] | 135 | 'e3': [3], 'e4':[4], 'e5': [5], |
---|
| 136 | 'e6': [6], 'e7': [7], 'e8': [8], |
---|
| 137 | 'e9': [9], 'e10': [10], 'e11': [11], |
---|
| 138 | 'e12': [12], 'e13': [13], 'e14': [14], |
---|
| 139 | 'e15': [15]}, |
---|
[3695] | 140 | 'maximum_triangle_area': 200000, |
---|
[3683] | 141 | 'filename': meshname, |
---|
[3679] | 142 | 'interior_regions': interior_regions}, |
---|
| 143 | verbose = True, evaluate=False) |
---|
[3559] | 144 | |
---|
| 145 | |
---|
| 146 | #------------------------------------------------------------------------------- |
---|
| 147 | # Setup computational domain |
---|
| 148 | #------------------------------------------------------------------------------- |
---|
[3679] | 149 | domain = Domain(meshname, use_cache = True, verbose = True) |
---|
[3559] | 150 | |
---|
| 151 | print 'Number of triangles = ', len(domain) |
---|
| 152 | print 'The extent is ', domain.get_extent() |
---|
| 153 | print domain.statistics() |
---|
| 154 | |
---|
| 155 | domain.set_name(project.basename) |
---|
| 156 | domain.set_datadir(project.outputtimedir) |
---|
| 157 | domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) |
---|
[3642] | 158 | domain.set_minimum_storable_height(0.01) |
---|
[3683] | 159 | domain.set_store_vertices_uniquely(False) # for writting to sww |
---|
[3559] | 160 | |
---|
| 161 | #------------------------------------------------------------------------------- |
---|
| 162 | # Setup initial conditions |
---|
| 163 | #------------------------------------------------------------------------------- |
---|
| 164 | |
---|
| 165 | tide = 0.0 |
---|
| 166 | |
---|
| 167 | domain.set_quantity('stage', tide) |
---|
| 168 | domain.set_quantity('friction', 0.0) |
---|
| 169 | |
---|
| 170 | domain.set_quantity('elevation', |
---|
| 171 | # filename = project.onshore_dem_name + '.pts', |
---|
[3683] | 172 | filename = project.combined_dem_name_2 + '.pts', |
---|
[3559] | 173 | # filename = project.offshore_dem_name + '.pts', |
---|
| 174 | use_cache = True, |
---|
| 175 | verbose = True, |
---|
| 176 | alpha = 0.1 |
---|
| 177 | ) |
---|
| 178 | |
---|
| 179 | #------------------------------------------------------------------------------- |
---|
[3661] | 180 | # Setup boundary conditions |
---|
[3559] | 181 | #------------------------------------------------------------------------------- |
---|
| 182 | print 'start ferret2sww' |
---|
[3626] | 183 | from anuga.shallow_water.data_manager import ferret2sww |
---|
[3679] | 184 | |
---|
[3559] | 185 | south = project.south |
---|
| 186 | north = project.north |
---|
| 187 | west = project.west |
---|
| 188 | east = project.east |
---|
| 189 | |
---|
| 190 | #note only need to do when an SWW file for the MOST boundary doesn't exist |
---|
| 191 | cache(ferret2sww, |
---|
| 192 | (source_dir + project.boundary_basename, |
---|
| 193 | source_dir + project.boundary_basename), |
---|
| 194 | {'verbose': True, |
---|
| 195 | 'minlat': south, |
---|
| 196 | 'maxlat': north, |
---|
| 197 | 'minlon': west, |
---|
| 198 | 'maxlon': east, |
---|
| 199 | # 'origin': project.mesh_origin, |
---|
| 200 | 'origin': domain.geo_reference.get_origin(), |
---|
| 201 | 'mean_stage': tide, |
---|
| 202 | 'zscale': 1, #Enhance tsunami |
---|
| 203 | 'fail_on_NaN': False, |
---|
| 204 | 'inverted_bathymetry': True}, |
---|
| 205 | #evaluate = True, |
---|
| 206 | verbose = True, |
---|
| 207 | dependencies = source_dir + project.boundary_basename + '.sww') |
---|
| 208 | |
---|
[3679] | 209 | |
---|
[3559] | 210 | print 'Available boundary tags', domain.get_boundary_tags() |
---|
| 211 | |
---|
[3679] | 212 | Bf = File_boundary(source_dir + project.boundary_basename + '.sww', |
---|
| 213 | domain, verbose = True) |
---|
[3559] | 214 | Br = Reflective_boundary(domain) |
---|
| 215 | Bd = Dirichlet_boundary([tide,0,0]) |
---|
| 216 | |
---|
| 217 | |
---|
| 218 | # 7 min square wave starting at 1 min, 6m high |
---|
| 219 | Bw = Time_boundary(domain = domain, |
---|
[3671] | 220 | f=lambda t: [(60<t<480)*10, 0, 0]) |
---|
[3559] | 221 | |
---|
[3615] | 222 | # for MOST BC |
---|
| 223 | #domain.set_boundary( {'top': Bd, 'left': Bd, |
---|
| 224 | # 'bottom': Bf, 'right': Bf} ) |
---|
[3559] | 225 | |
---|
[3615] | 226 | # for testing |
---|
[3673] | 227 | #domain.set_boundary( {'topr': Bd, 'left': Bd, 'top': Bd, |
---|
| 228 | # 'bottom': Bw, 'bright': Bd} ) |
---|
| 229 | domain.set_boundary( {'e0': Bd, 'e1': Bd, 'e2': Bd, 'e3': Bd, 'e4': Bd, |
---|
| 230 | 'e5': Bd, 'e6': Bd, 'e7': Bd, 'e8': Bd, 'e9': Bd, |
---|
[3683] | 231 | 'e10': Bd, 'e11': Bd, 'e12': Bf, 'e13': Bf, 'e14': Bf, |
---|
| 232 | 'e15': Bf} ) |
---|
[3615] | 233 | |
---|
[3559] | 234 | #------------------------------------------------------------------------------- |
---|
| 235 | # Evolve system through time |
---|
| 236 | #------------------------------------------------------------------------------- |
---|
| 237 | import time |
---|
| 238 | t0 = time.time() |
---|
| 239 | |
---|
[3679] | 240 | for t in domain.evolve(yieldstep = 240, finaltime = 6800): |
---|
[3559] | 241 | domain.write_time() |
---|
[3615] | 242 | domain.write_boundary_statistics(tags = 'bottom') |
---|
[3559] | 243 | |
---|
[3679] | 244 | for t in domain.evolve(yieldstep = 30, finaltime = 9000 |
---|
| 245 | ,skip_initial_step = True): |
---|
| 246 | domain.write_time() |
---|
| 247 | domain.write_boundary_statistics(tags = 'bottom') |
---|
| 248 | |
---|
| 249 | for t in domain.evolve(yieldstep = 240, finaltime = 15000 |
---|
| 250 | ,skip_initial_step = True): |
---|
| 251 | domain.write_time() |
---|
| 252 | domain.write_boundary_statistics(tags = 'bottom') |
---|
[3559] | 253 | |
---|
| 254 | print 'That took %.2f seconds' %(time.time()-t0) |
---|
| 255 | |
---|
| 256 | print 'finished' |
---|