[3477] | 1 | """Script for running a tsunami inundation scenario for Broome, WA, 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, |
---|
| 8 | the elevation data and a simulated submarine landslide. |
---|
| 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 | from os import sep |
---|
| 17 | from os.path import dirname, basename |
---|
| 18 | import time |
---|
| 19 | |
---|
| 20 | # Related major packages |
---|
[3514] | 21 | from anuga.pyvolution.shallow_water import Domain, Reflective_boundary, \ |
---|
[3477] | 22 | Dirichlet_boundary, Time_boundary, File_boundary |
---|
[3514] | 23 | from anuga.pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts |
---|
| 24 | from anuga.pyvolution.combine_pts import combine_rectangular_points_files |
---|
| 25 | from anuga.pyvolution.pmesh2domain import pmesh_to_domain_instance |
---|
[3477] | 26 | from shutil import copy |
---|
| 27 | from os import mkdir, access, F_OK |
---|
[3514] | 28 | from anuga.geospatial_data.geospatial_data import * |
---|
[3477] | 29 | import sys |
---|
[3514] | 30 | from anuga.pyvolution.util import Screen_Catcher |
---|
[3477] | 31 | |
---|
| 32 | # Application specific imports |
---|
| 33 | import 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 if dir doesn't exist |
---|
| 41 | if access(project.outputtimedir,F_OK) == 0 : |
---|
| 42 | mkdir (project.outputtimedir) |
---|
| 43 | copy (dirname(project.__file__) +sep+ project.__name__+'.py', project.outputtimedir + project.__name__+'.py') |
---|
| 44 | copy (__file__, project.outputtimedir + basename(__file__)) |
---|
| 45 | print 'project.outputtimedir',project.outputtimedir |
---|
| 46 | |
---|
| 47 | # normal screen output is stored in |
---|
| 48 | screen_output_name = project.outputtimedir + "screen_output.txt" |
---|
| 49 | screen_error_name = project.outputtimedir + "screen_error.txt" |
---|
| 50 | |
---|
| 51 | # used to catch screen output to file |
---|
| 52 | sys.stdout = Screen_Catcher(screen_output_name) |
---|
| 53 | sys.stderr = Screen_Catcher(screen_error_name) |
---|
| 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 | # Do for coarse and fine data |
---|
| 61 | # Fine pts file to be clipped to area of interest |
---|
| 62 | #------------------------------------------------------------------------------- |
---|
| 63 | |
---|
| 64 | # filenames |
---|
| 65 | meshname = project.meshname+'.msh' |
---|
| 66 | source_dir = project.boundarydir |
---|
| 67 | |
---|
| 68 | # creates DEM from asc data |
---|
| 69 | convert_dem_from_ascii2netcdf(project.onshore_dem_name, use_cache=True, verbose=True) |
---|
| 70 | |
---|
| 71 | #creates pts file from DEM |
---|
| 72 | dem2pts(project.onshore_dem_name, |
---|
| 73 | easting_min=project.eastingmin, |
---|
| 74 | easting_max=project.eastingmax, |
---|
| 75 | northing_min=project.northingmin, |
---|
| 76 | northing_max= project.northingmax, |
---|
| 77 | use_cache=True, |
---|
| 78 | verbose=True) |
---|
| 79 | |
---|
| 80 | print 'create G1' |
---|
| 81 | G1 = Geospatial_data(file_name = project.offshore_dem_name1 + '.xya') |
---|
| 82 | print 'create G2' |
---|
| 83 | G2 = Geospatial_data(file_name = project.offshore_dem_name2 + '.xya') |
---|
| 84 | print 'create G3' |
---|
| 85 | G3 = Geospatial_data(file_name = project.onshore_dem_name + '.pts') |
---|
| 86 | print 'create G4' |
---|
| 87 | G4 = Geospatial_data(file_name = project.coast_dem_name + '.xya') |
---|
| 88 | print 'add G1+G2+G3+G4' |
---|
| 89 | G = G1 + G2 + G3 + G4 |
---|
| 90 | print 'export G' |
---|
| 91 | G.export_points_file(project.combined_dem_name + '.pts') |
---|
| 92 | |
---|
| 93 | #------------------------------------------------------------------------------- |
---|
| 94 | # Create the triangular mesh based on overall clipping polygon with a tagged |
---|
| 95 | # boundary and interior regions defined in project.py along with |
---|
| 96 | # resolutions (maximal area of per triangle) for each polygon |
---|
| 97 | #------------------------------------------------------------------------------- |
---|
| 98 | |
---|
| 99 | from pmesh.mesh_interface import create_mesh_from_regions |
---|
| 100 | |
---|
| 101 | region_res = 500000 |
---|
| 102 | coast_res = 500 |
---|
| 103 | broome_res = 5000 |
---|
| 104 | interior_regions = [[project.poly_broome, broome_res], |
---|
| 105 | [project.poly_region, region_res]] |
---|
| 106 | |
---|
| 107 | print 'number of interior regions', len(interior_regions) |
---|
| 108 | |
---|
[3514] | 109 | from anuga.utilities.polygon import plot_polygons |
---|
[3477] | 110 | if sys.platform == 'win32': |
---|
| 111 | #figname = project.outputtimedir + 'pt_hedland_polys' |
---|
| 112 | figname = 'broome_polys_test' |
---|
| 113 | plot_polygons([project.polyAll,project.poly_broome,project.poly_region], |
---|
| 114 | figname, |
---|
| 115 | verbose = True) |
---|
| 116 | |
---|
| 117 | print 'start create mesh from regions' |
---|
| 118 | from caching import cache |
---|
| 119 | _ = cache(create_mesh_from_regions, |
---|
| 120 | project.polyAll, |
---|
| 121 | {'boundary_tags': {'topright': [0], 'topleft': [1], |
---|
| 122 | 'left': [2], 'bottom0': [3], |
---|
| 123 | 'bottom1': [4], 'bottom2': [5], |
---|
| 124 | 'bottom3': [6], 'right': [7]}, |
---|
| 125 | 'maximum_triangle_area': 250000, |
---|
| 126 | 'filename': meshname, |
---|
| 127 | 'interior_regions': interior_regions}, |
---|
| 128 | verbose = True, evaluate=True) |
---|
| 129 | |
---|
| 130 | #------------------------------------------------------------------------------- |
---|
| 131 | # Setup computational domain |
---|
| 132 | #------------------------------------------------------------------------------- |
---|
| 133 | domain = Domain(meshname, use_cache = False, verbose = True) |
---|
| 134 | |
---|
| 135 | print domain.statistics() |
---|
| 136 | print 'Number of triangles = ', len(domain) |
---|
| 137 | print 'The extent is ', domain.get_extent() |
---|
| 138 | print domain.statistics() |
---|
| 139 | |
---|
| 140 | domain.set_name(project.basename) |
---|
| 141 | domain.set_datadir(project.outputtimedir) |
---|
| 142 | domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) |
---|
| 143 | |
---|
| 144 | #------------------------------------------------------------------------------- |
---|
| 145 | # Setup initial conditions |
---|
| 146 | #------------------------------------------------------------------------------- |
---|
| 147 | |
---|
| 148 | tide = 0. |
---|
| 149 | #high |
---|
| 150 | #tide = |
---|
| 151 | #low |
---|
| 152 | #tide = |
---|
| 153 | |
---|
| 154 | domain.set_quantity('stage', tide) |
---|
| 155 | domain.set_quantity('friction', 0.0) |
---|
| 156 | print 'hi and file',project.combined_dem_name + '.pts' |
---|
| 157 | |
---|
| 158 | domain.set_quantity('elevation', |
---|
| 159 | filename = project.combined_dem_name + '.pts', |
---|
| 160 | use_cache = True, |
---|
| 161 | verbose = True, |
---|
| 162 | alpha = 0.1 |
---|
| 163 | ) |
---|
| 164 | |
---|
| 165 | #------------------------------------------------------------------------------- |
---|
| 166 | # Setup boundary conditions |
---|
| 167 | #------------------------------------------------------------------------------- |
---|
| 168 | print 'start ferret2sww' |
---|
| 169 | # skipped as results in file SU-AU_clipped is correct for all WA |
---|
| 170 | |
---|
[3514] | 171 | from anuga.pyvolution.data_manager import ferret2sww |
---|
[3477] | 172 | |
---|
| 173 | south = project.south |
---|
| 174 | north = project.north |
---|
| 175 | west = project.west |
---|
| 176 | east = project.east |
---|
| 177 | |
---|
| 178 | #note only need to do when an SWW file for the MOST boundary doesn't exist |
---|
| 179 | cache(ferret2sww, |
---|
| 180 | (source_dir + project.boundary_basename, |
---|
| 181 | source_dir + project.boundary_basename+'_'+project.basename), |
---|
| 182 | {'verbose': True, |
---|
| 183 | 'minlat': south, |
---|
| 184 | 'maxlat': north, |
---|
| 185 | 'minlon': west, |
---|
| 186 | 'maxlon': east, |
---|
| 187 | # 'origin': project.mesh_origin, |
---|
| 188 | 'origin': domain.geo_reference.get_origin(), |
---|
| 189 | 'mean_stage': tide, |
---|
| 190 | 'zscale': 1, #Enhance tsunami |
---|
| 191 | 'fail_on_NaN': False, |
---|
| 192 | 'inverted_bathymetry': True}, |
---|
| 193 | evaluate = True, |
---|
| 194 | verbose = True, |
---|
| 195 | dependencies = source_dir + project.boundary_basename + '.sww') |
---|
| 196 | |
---|
| 197 | print 'Available boundary tags', domain.get_boundary_tags() |
---|
| 198 | |
---|
| 199 | Bf = File_boundary(source_dir + project.boundary_basename + '.sww', |
---|
| 200 | domain, verbose = True) |
---|
| 201 | Br = Reflective_boundary(domain) |
---|
| 202 | Bd = Dirichlet_boundary([tide,0,0]) |
---|
| 203 | domain.set_boundary( {'topright': Bf,'topleft': Bf, 'left': Bd, 'bottom0': Bd, |
---|
| 204 | 'bottom1': Bd, 'bottom2': Bd, 'bottom3': Bd, |
---|
| 205 | 'right': Bd}) |
---|
| 206 | |
---|
| 207 | #------------------------------------------------------------------------------- |
---|
| 208 | # Evolve system through time |
---|
| 209 | #------------------------------------------------------------------------------- |
---|
| 210 | import time |
---|
| 211 | t0 = time.time() |
---|
| 212 | |
---|
| 213 | for t in domain.evolve(yieldstep = 240, finaltime = 10800): |
---|
| 214 | domain.write_time() |
---|
| 215 | domain.write_boundary_statistics(tags = 'topright') |
---|
| 216 | |
---|
| 217 | for t in domain.evolve(yieldstep = 120, finaltime = 16200 |
---|
| 218 | ,skip_initial_step = True): |
---|
| 219 | domain.write_time() |
---|
| 220 | domain.write_boundary_statistics(tags = 'topright') |
---|
| 221 | |
---|
| 222 | for t in domain.evolve(yieldstep = 60, finaltime = 21600 |
---|
| 223 | ,skip_initial_step = True): |
---|
| 224 | domain.write_time() |
---|
| 225 | domain.write_boundary_statistics(tags = 'topright') |
---|
| 226 | |
---|
| 227 | for t in domain.evolve(yieldstep = 120, finaltime = 27000 |
---|
| 228 | ,skip_initial_step = True): |
---|
| 229 | domain.write_time() |
---|
| 230 | domain.write_boundary_statistics(tags = 'topright') |
---|
| 231 | |
---|
| 232 | for t in domain.evolve(yieldstep = 240, finaltime = 36000 |
---|
| 233 | ,skip_initial_step = True): |
---|
| 234 | domain.write_time() |
---|
| 235 | domain.write_boundary_statistics(tags = 'topright') |
---|
| 236 | |
---|
| 237 | print 'That took %.2f seconds' %(time.time()-t0) |
---|
| 238 | |
---|
| 239 | print 'finished' |
---|