[2461] | 1 | """Script for running a tsunami inundation scenario for Sydney, NSW, 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 a simulated submarine landslide. |
---|
| 9 | |
---|
| 10 | """ |
---|
| 11 | |
---|
| 12 | |
---|
[2887] | 13 | #------------------------------------------------------------------------------ |
---|
[2478] | 14 | # Import necessary modules |
---|
[2887] | 15 | #------------------------------------------------------------------------------ |
---|
[2461] | 16 | |
---|
| 17 | # Standard modules |
---|
| 18 | import os |
---|
| 19 | import time |
---|
| 20 | |
---|
| 21 | # Related major packages |
---|
[3563] | 22 | from anuga.shallow_water import Domain, Dirichlet_boundary |
---|
| 23 | from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts |
---|
| 24 | from anuga.abstract_2d_finite_volumes.combine_pts import combine_rectangular_points_files |
---|
[3534] | 25 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
[2461] | 26 | |
---|
| 27 | # Application specific imports |
---|
[2887] | 28 | import project # Define file names and polygons |
---|
[3563] | 29 | from anuga.shallow_water.smf import slump_tsunami # Function for submarine mudslide |
---|
[2461] | 30 | |
---|
| 31 | |
---|
| 32 | |
---|
[2887] | 33 | #------------------------------------------------------------------------------ |
---|
[2599] | 34 | # Prepare bathymetric and topographic data |
---|
[2887] | 35 | #------------------------------------------------------------------------------ |
---|
[2478] | 36 | |
---|
[2461] | 37 | # filenames |
---|
[3869] | 38 | #coarsedemname = project.coarsedemname + '.pts' |
---|
| 39 | #finedemname = project.finedemname + '.pts' |
---|
| 40 | demname = project.demname |
---|
[2461] | 41 | meshname = project.meshname+'.msh' |
---|
| 42 | |
---|
[3869] | 43 | convert_dem_from_ascii2netcdf(demname, use_cache=True, verbose=True) |
---|
[2461] | 44 | |
---|
[3869] | 45 | # creates pts file |
---|
| 46 | dem2pts(demname, use_cache=True, verbose=True) |
---|
| 47 | |
---|
[2461] | 48 | # combining the coarse and fine data |
---|
[3869] | 49 | #combine_rectangular_points_files(finedemname, |
---|
| 50 | # coarsedemname, |
---|
| 51 | # combineddemname) |
---|
[2461] | 52 | |
---|
| 53 | |
---|
[2478] | 54 | |
---|
[2887] | 55 | #------------------------------------------------------------------------------ |
---|
[2478] | 56 | # Setup computational domain |
---|
[2887] | 57 | #------------------------------------------------------------------------------ |
---|
[2478] | 58 | |
---|
| 59 | # Interior regions and mesh resolutions |
---|
[2461] | 60 | interior_res = 5000 |
---|
[3944] | 61 | shallow_res = 15000 |
---|
| 62 | ##interior_regions = [[project.northern_polygon, interior_res], |
---|
| 63 | ## [project.harbour_polygon, interior_res], |
---|
| 64 | ## [project.southern_polygon, interior_res], |
---|
| 65 | ## [project.manly_polygon, interior_res], |
---|
| 66 | ## [project.top_polygon, interior_res]] |
---|
| 67 | interior_regions = [[project.coastal_polygon, interior_res], |
---|
| 68 | [project.shallow_polygon, shallow_res]] |
---|
[2461] | 69 | |
---|
[3136] | 70 | create_mesh_from_regions(project.demopoly, |
---|
| 71 | boundary_tags= {'oceannorth': [0], |
---|
[3190] | 72 | 'onshorenorth1': [1], |
---|
| 73 | 'onshorenorthwest1': [2], |
---|
| 74 | 'onshorenorthwest2': [3], |
---|
| 75 | 'onshorenorth2': [4], |
---|
| 76 | 'onshorewest': [5], |
---|
| 77 | 'onshoresouth': [6], |
---|
| 78 | 'oceansouth': [7], |
---|
| 79 | 'oceaneast': [8]}, |
---|
[2478] | 80 | maximum_triangle_area=100000, |
---|
| 81 | filename=meshname, |
---|
| 82 | interior_regions=interior_regions) |
---|
[2461] | 83 | |
---|
| 84 | |
---|
[2805] | 85 | #Create shallow water domain |
---|
[2461] | 86 | |
---|
[2866] | 87 | domain = Domain(meshname, |
---|
| 88 | use_cache = True, |
---|
| 89 | verbose = True) |
---|
[2805] | 90 | |
---|
[2866] | 91 | |
---|
[2461] | 92 | print 'Number of triangles = ', len(domain) |
---|
| 93 | print 'The extent is ', domain.get_extent() |
---|
| 94 | |
---|
| 95 | domain.set_name(project.basename) |
---|
| 96 | domain.set_datadir(project.outputdir) |
---|
| 97 | domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) |
---|
| 98 | |
---|
| 99 | |
---|
[2887] | 100 | #------------------------------------------------------------------------------ |
---|
[2461] | 101 | # Set up scenario (tsunami_source is a callable object used with set_quantity) |
---|
[2887] | 102 | #------------------------------------------------------------------------------ |
---|
[2461] | 103 | |
---|
| 104 | tsunami_source = slump_tsunami(length=30000.0, |
---|
| 105 | depth=400.0, |
---|
| 106 | slope=6.0, |
---|
| 107 | thickness=176.0, |
---|
| 108 | radius=3330, |
---|
| 109 | dphi=0.23, |
---|
| 110 | x0=project.slump_origin[0], |
---|
| 111 | y0=project.slump_origin[1], |
---|
| 112 | alpha=0.0, |
---|
| 113 | domain=domain) |
---|
| 114 | |
---|
| 115 | |
---|
[2887] | 116 | #------------------------------------------------------------------------------ |
---|
[2461] | 117 | # Setup initial conditions |
---|
[2887] | 118 | #------------------------------------------------------------------------------ |
---|
[2461] | 119 | |
---|
| 120 | domain.set_quantity('stage', tsunami_source) |
---|
| 121 | domain.set_quantity('friction', 0.0) |
---|
| 122 | domain.set_quantity('elevation', |
---|
[3869] | 123 | filename = demname + '.pts', |
---|
[2461] | 124 | use_cache = True, |
---|
| 125 | verbose = True) |
---|
| 126 | |
---|
[2887] | 127 | #------------------------------------------------------------------------------ |
---|
[3136] | 128 | # Setup boundary conditions (all Dirichlet) |
---|
[2887] | 129 | #------------------------------------------------------------------------------ |
---|
[2461] | 130 | |
---|
| 131 | print 'Available boundary tags', domain.get_boundary_tags() |
---|
| 132 | |
---|
[3136] | 133 | Bd = Dirichlet_boundary([0.0,0.0,0.0]) |
---|
[3190] | 134 | domain.set_boundary( {'oceannorth': Bd, |
---|
| 135 | 'onshorenorth1': Bd, |
---|
| 136 | 'onshorenorthwest1': Bd, |
---|
| 137 | 'onshorenorthwest2': Bd, |
---|
| 138 | 'onshorenorth2': Bd, |
---|
| 139 | 'onshorewest': Bd, |
---|
| 140 | 'onshoresouth': Bd, |
---|
| 141 | 'oceansouth': Bd, |
---|
| 142 | 'oceaneast': Bd} ) |
---|
[2461] | 143 | |
---|
[2887] | 144 | #------------------------------------------------------------------------------ |
---|
[2461] | 145 | # Evolve system through time |
---|
[2887] | 146 | #------------------------------------------------------------------------------ |
---|
[2461] | 147 | |
---|
| 148 | import time |
---|
| 149 | t0 = time.time() |
---|
| 150 | |
---|
[3190] | 151 | for t in domain.evolve(yieldstep = 60, finaltime = 7200): |
---|
[2497] | 152 | print domain.timestepping_statistics() |
---|
[3136] | 153 | print domain.boundary_statistics(tags = 'oceaneast') |
---|
[2461] | 154 | |
---|
| 155 | print 'That took %.2f seconds' %(time.time()-t0) |
---|