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