Changeset 5194
- Timestamp:
- Apr 9, 2008, 9:09:18 AM (17 years ago)
- Location:
- anuga_work/production/tonga
- Files:
-
- 7 added
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/production/tonga/project.py
r5139 r5194 30 30 #Making assumptions about the location of scenario data 31 31 state = 'sw_pacific' 32 scenario_name = 'Tonga '33 scenario = ' fixed_wave'32 scenario_name = 'Tonga_slide' 33 scenario = 'tonga' 34 34 35 35 tide = 0 … … 42 42 export_cellsize=50 43 43 setup='trial' 44 source='test' 44 45 45 46 if setup =='trial': … … 59 60 yieldstep=60 60 61 61 dir_comment='_'+setup+'_'+str(tide)+'_'+str(source)+'_'+str(user) 62 63 64 topo_gridfile ='tongatapu_10mgrid' 62 dir_comment='_'+setup+'_'+str(tide)+'_'+str(scenario_name)+'_'+str(user) 63 64 # onshore data 5m countour 65 onshore_name = 'topography_tongatapu' # original' 66 67 # AHO + DPI data + colin French coastline 68 #coast_name = 'waterline' 69 Singlebeam_name = 'Tongatapu_SB_5m grid' 70 Multibeam_name = 'Tongatapu_MB_30m grid' 71 Chart_name= 'Tongatapu_Chart' 72 Derived_bath_name= 'Derived_Bathy' 73 added_data_name='joining_eastIsland_toreef' 74 75 76 #final topo name 77 combined_name ='Tongatapu_combined_elevation' 78 combined_smaller_name = 'Tongatapu_combined_elevation_smaller' 65 79 66 80 anuga_dir = home+state+sep+scenario+sep+'anuga'+sep … … 70 84 #topographies_time_dir = topographies_dir+build_time+sep 71 85 86 # input topo file location 87 onshore_in_dir_name = topographies_in_dir + onshore_name 88 89 Singlebeam_in_dir_name = topographies_in_dir + Singlebeam_name 90 Multibeam_in_dir_name = topographies_in_dir + Multibeam_name 91 Chart_in_dir_name = topographies_in_dir + Chart_name 92 Derived_bath_in_dir_name = topographies_in_dir + Derived_bath_name 93 added_data_in_dir_name = topographies_in_dir + added_data_name 94 95 onshore_dir_name = topographies_dir + onshore_name 96 97 Singlebeam_dir_name = topographies_dir + Singlebeam_name 98 Multibeam_dir_name = topographies_dir + Multibeam_name 99 Chart_dir_name = topographies_dir + Chart_name 100 Derived_bath_dir_name = topographies_dir + Derived_bath_name 101 added_data_dir_name = topographies_dir + added_data_name 102 103 #final topo files 104 combined_dir_name = topographies_dir + combined_name 105 #combined_time_dir_name = topographies_time_dir + combined_name 106 combined_smaller_name_dir = topographies_dir + combined_smaller_name 107 #combined_time_dir_final_name = topographies_time_dir + combined_final_name 108 72 109 meshes_dir = anuga_dir+'meshes'+sep 73 110 meshes_dir_name = meshes_dir + scenario_name … … 76 113 tide_dir = anuga_dir+'tide_data'+sep 77 114 115 116 #boundaries_source = '1' 117 118 if source =='dampier': 119 boundaries_name = 'broome_3854_17042007' #Dampier gun 120 boundaries_in_dir = anuga_dir+'boundaries'+sep+sep+'dampier'+sep+'1_10000'+sep 121 122 if source=='onslow': 123 boundaries_name = 'broome_3859_16052007' #onslow_hedland_broome gun 124 boundaries_in_dir = anuga_dir+'boundaries'+sep+sep+'onslow_hedland_broome'+sep+'1_10000'+sep 125 126 if source=='exmouth': 127 boundaries_name = 'broome_3103_18052007' #exmouth gun 128 boundaries_in_dir = anuga_dir+'boundaries'+sep+sep+'exmouth'+sep+'1_10000'+sep 129 130 if source=='test': 131 boundaries_name = 'other' #exmouth gun 132 boundaries_in_dir = anuga_dir+'boundaries'+sep 133 134 135 #boundaries locations 78 136 boundaries_in_dir_name = boundaries_in_dir + scenario_name 79 137 boundaries_dir = anuga_dir+'boundaries'+sep 80 138 boundaries_dir_name = boundaries_dir + scenario_name 139 #boundaries_time_dir = anuga_dir+'boundaries'+sep+build_time+sep 140 #boundaries_time_dir_name = boundaries_time_dir + boundaries_name #Used by post processing 81 141 82 142 #output locations … … 84 144 output_build_time_dir = anuga_dir+'outputs'+sep+build_time+dir_comment+sep 85 145 output_run_time_dir = anuga_dir+'outputs'+sep+run_time+dir_comment+sep 86 output_run_time_dir_name = output_run_time_dir + scenario_name #Used by post pr 146 output_run_time_dir_name = output_run_time_dir + scenario_name #Used by post processing 147 148 #gauges 149 #gauge_name = 'perth.csv' 150 #gauges_dir = anuga_dir+'gauges'+sep 151 #gauges_dir_name = gauges_dir + gauge_name 152 153 #buildings_filename = gauges_dir + 'Perth_res_Project.csv' 154 #buildings_filename_out = 'Perth_res_Project_modified.csv' 87 155 88 156 ############################### … … 92 160 93 161 poly_all = read_polygon(polygons_dir+'extent.txt') 94 res_poly_all = 100000*res_factor 95 162 res_poly_all = 40000*res_factor 163 164 #refzone = 50 96 165 97 166 ############################### … … 99 168 ############################### 100 169 101 # interior polygons 102 poly_tongatapu = read_polygon('tongatapu_finemesharea.txt') 103 poly_island1 = read_polygon('island1.txt') 104 poly_island2 = read_polygon('island2.txt') 105 poly_island3 = read_polygon('islands3.txt') 106 poly_island4 = read_polygon('islands4.txt') 107 poly_island5= read_polygon('islands5.txt') 108 poly_island6= read_polygon('islands6.txt') 109 poly_island7= read_polygon('islands7.txt') 110 poly_island8= read_polygon('islands8.txt') 111 poly_island9= read_polygon('islands9.txt') 112 poly_island10= read_polygon('islands10.txt') 113 poly_island11= read_polygon('islands11.txt') 114 poly_island12= read_polygon('islands12.txt') 115 poly_island13= read_polygon('islands13.txt') 170 poly_tongatapu = read_polygon(polygons_dir+'finres.txt') 171 poly_island1 = read_polygon(polygons_dir+'island1.txt') 172 poly_island2 = read_polygon(polygons_dir+'island2.txt') 173 poly_island3 = read_polygon(polygons_dir+'island3.txt') 174 poly_island4 = read_polygon(polygons_dir+'island4.txt') 175 poly_island5= read_polygon(polygons_dir+'island5.txt') 176 poly_island6= read_polygon(polygons_dir+'island6.txt') 177 poly_island7= read_polygon(polygons_dir+'island7.txt') 178 poly_island8= read_polygon(polygons_dir+'island8.txt') 179 poly_island9= read_polygon(polygons_dir+'island9.txt') 180 poly_island10= read_polygon(polygons_dir+'island10.txt') 181 poly_island11= read_polygon(polygons_dir+'island11.txt') 182 poly_island12= read_polygon(polygons_dir+'island12.txt') 183 poly_island13= read_polygon(polygons_dir+'island13.txt') 116 184 117 185 118 186 islands_res = 5000 119 cairns_res = 5000 120 interior_regions = [[project.poly_tongatapu, tongatapu_res], 121 [project.poly_island1, islands_res], 122 [project.poly_island2, islands_res], 123 [project.poly_island3, islands_res], 124 [project.poly_island4, islands_res], 125 [project.poly_island5, islands_res], 126 [project.poly_island6, islands_res], 127 [project.poly_island7, islands_res], 128 [project.poly_island8, islands_res], 129 [project.poly_island9, islands_res], 130 [project.poly_island10, islands_res], 131 [project.poly_island11, islands_res], 132 [project.poly_island12, islands_res], 133 [project.poly_island13, islands_res]] 134 135 boundary_tags={'ocean_east': [0],'ocean_north': [1],'ocean_west': [2],'land': [3,4,5,6,7],'ocean_southeast':[8,9]} 187 tongatapu_res = 5000 188 interior_regions = [[poly_tongatapu, tongatapu_res], 189 [poly_island1, islands_res], 190 [poly_island2, islands_res], 191 [poly_island3, islands_res], 192 [poly_island4, islands_res], 193 [poly_island5, islands_res], 194 [poly_island6, islands_res], 195 [poly_island7, islands_res], 196 [poly_island8, islands_res], 197 [poly_island9, islands_res], 198 [poly_island10, islands_res], 199 [poly_island11, islands_res], 200 [poly_island12, islands_res], 201 [poly_island13, islands_res]] 202 203 boundary_tags={'ocean_east':[0], 204 'ocean_north':[1], 205 'ocean_west':[2,3], 206 'land':[4,5,6,7,8,9,10,11,12,13,14], 207 'ocean_southeast':[15,16]} 136 208 137 209 trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all) 138 210 139 211 print 'min number triangles', trigs_min 140 141 212 ################################################################### 142 213 # Clipping regions for export to asc and regions for clipping data … … 156 227 157 228 158 159 160 161 162 163 164 165 # topo Filenames166 dem_name = 'tongatapu_10mgrid'167 meshname = 'tongatapu.msh'168 169 # Create DEM from asc data170 convert_dem_from_ascii2netcdf(dem_name, use_cache=True, verbose=True)171 172 # Create pts file for onshore DEM173 dem2pts(dem_name, use_cache=True, verbose=True)174 175 176 177 178 179 180 # -*- coding: cp1252 -*-181 """Common filenames and locations for topographic data, meshes and outputs.182 """183 184 from os import sep, environ, getenv, getcwd185 from os.path import expanduser186 import sys187 from time import localtime, strftime, gmtime188 from anuga.utilities.polygon import read_polygon, plot_polygons, \189 polygon_area, is_inside_polygon190 191 ###############################192 # Domain definitions193 ###############################194 195 # bounding polygon for study area196 bounding_polygon = read_polygon('Y:\data\sw_pacific\tonga\anuga\boundaries\extent.txt')197 198 print 'Area of bounding polygon', polygon_area(bounding_polygon)/1000000.0199 200 ###############################201 # Interior region definitions202 ###############################203 204 # interior polygons205 poly_tongatapu = read_polygon('tongatapu_finemesharea.txt')206 poly_island1 = read_polygon('island1.txt')207 poly_island2 = read_polygon('island2.txt')208 poly_island3 = read_polygon('islands3.txt')209 poly_island4 = read_polygon('islands4.txt')210 poly_island5= read_polygon('islands5.txt')211 poly_island6= read_polygon('islands6.txt')212 poly_island7= read_polygon('islands7.txt')213 poly_island8= read_polygon('islands8.txt')214 poly_island9= read_polygon('islands9.txt')215 poly_island10= read_polygon('islands10.txt')216 poly_island11= read_polygon('islands11.txt')217 poly_island12= read_polygon('islands12.txt')218 poly_island13= read_polygon('islands13.txt')219 220 #plot_polygons([bounding_polygon,poly_cairns,poly_island0,poly_island1,\221 # poly_island2,poly_island3,poly_shallow],\222 # 'boundingpoly',verbose=False)223 224 ###################################################################225 # Clipping regions for export to asc and regions for clipping data226 ###################################################################227 228 # exporting asc grid229 eastingmin = 670500230 eastingmax = 712750231 northingmin = 7646000232 northingmax = 7677000233 234 235 slide_origin = [701290, 7665750] # move onto the continental shelf, depth = 500236 slide_depth = 207.237 238 #gauge_filename = 'gauges.csv' -
anuga_work/production/tonga/tongatapu.py
r5139 r5194 24 24 25 25 # Related major packages 26 from anuga.shallow_water import Time_boundary 26 27 from anuga.shallow_water import Domain 27 28 from anuga.shallow_water import Dirichlet_boundary … … 31 32 from Numeric import allclose 32 33 from anuga.shallow_water.data_manager import export_grid 33 34 from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts 34 35 from anuga.pmesh.mesh_interface import create_mesh_from_regions 35 36 from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters … … 81 82 # resolutions (maximal area of per triangle) for each polygon 82 83 #-------------------------------------------------------------------------- 83 84 #IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it85 # causes problems with the ability to cache set quantity which takes alot of times86 84 87 85 if myid == 0: … … 98 96 barrier() 99 97 98 scenario='fixed_wave' 99 100 100 #------------------------------------------------------------------------- 101 101 # Setup computational domain … … 133 133 domain.set_quantity('elevation', 134 134 filename = kwargs['bathy_file'], 135 use_cache = True,135 use_cache = False, 136 136 verbose = True, 137 137 alpha = kwargs['alpha']) … … 170 170 tsunami_source = slide_tsunami(length=35000.0, 171 171 depth=project.slide_depth, 172 slope= 6.0,173 thickness= 500.0,172 slope=10.0, 173 thickness=800.0, 174 174 x0=project.slide_origin[0], 175 175 y0=project.slide_origin[1], … … 192 192 # FIXME (Ole): Change this to Transmissive Momentum as 193 193 # suggested by Rajaraman (ticket:208) 194 from math import sin, pi 194 195 Bw = Time_boundary(domain = domain, 195 f=lambda t: [ (60<t<3660)*100, 0, 0])196 domain.set_boundary({'ocean_east': 197 'ocean_ southeast':Bd,198 ' land':Bd,199 ' ocean_north':Bd,200 'ocean_ west':})196 f=lambda t: [5.0*sin(t*pi*1/15), 0, 0]) 197 domain.set_boundary({'ocean_east':Bw, 198 'ocean_north':Bd, 199 'ocean_west':Bd, 200 'land':Bd, 201 'ocean_southeast':Bd}) 201 202 202 203 # boundary conditions for slide scenario 203 204 if scenario == 'slide': 204 domain.set_boundary({'ocean_east': 205 'ocean_ southeast':Bd,206 ' land':Bd,207 ' ocean_north':Bd,208 'ocean_ west':})205 domain.set_boundary({'ocean_east':Bd, 206 'ocean_north':Bd, 207 'ocean_west':Bd, 208 'land':Bd, 209 'ocean_southeast':Bd}) 209 210 210 211 #---------------------------------------------------------------------------- … … 213 214 t0 = time.time() 214 215 215 for t in domain.evolve(yieldstep = 240, finaltime = kwargs['finaltime']):216 for t in domain.evolve(yieldstep = 1, finaltime = kwargs['finaltime']): 216 217 domain.write_time() 217 218 domain.write_boundary_statistics(tags = 'ocean_east') … … 240 241 print 'memory usage after del domain',mem_usage() 241 242 242 swwfile = kwargs['output_dir']+kwargs[' scenario_name']243 swwfile = kwargs['output_dir']+kwargs['aa_scenario_name'] 243 244 print'swwfile',swwfile 244 245 245 export_grid(swwfile, extra_name_out = 'to wn',246 export_grid(swwfile, extra_name_out = 'tonga', 246 247 quantities = ['speed','depth','elevation','stage'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation 247 248 #quantities = ['speed','depth'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation 248 249 timestep = None, 249 250 reduction = max, 250 cellsize = kwargs['export_cellsize'],251 cellsize = 25, 251 252 NODATA_value = -1E-030, 252 253 easting_min = project.eastingmin, … … 271 272 kwargs['midtime']=project.midtime 272 273 kwargs['finaltime']=project.finaltime 273 274 274 kwargs['output_dir']=project.output_run_time_dir 275 275 kwargs['bathy_file']=project.combined_dir_name+'.txt'
Note: See TracChangeset
for help on using the changeset viewer.