[4832] | 1 | |
---|
| 2 | |
---|
| 3 | from anuga.geospatial_data.geospatial_data import * |
---|
[4846] | 4 | from os import sep, environ, getenv, getcwd,umask |
---|
[4832] | 5 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
| 6 | from anuga.shallow_water import Domain |
---|
| 7 | from anuga.shallow_water import Reflective_boundary |
---|
| 8 | from anuga.shallow_water.data_manager import export_grid |
---|
[4846] | 9 | from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files |
---|
[4832] | 10 | from math import sqrt |
---|
[4846] | 11 | from time import localtime, strftime, gmtime |
---|
[4832] | 12 | |
---|
[4846] | 13 | umask(002) |
---|
| 14 | time = strftime('%Y%m%d_%H%M%S',gmtime()) |
---|
| 15 | |
---|
| 16 | home = getenv('INUNDATIONHOME') + sep +'data'+sep |
---|
| 17 | output_dir = home + 'anuga_validation'+sep+'cocos'+sep+'anuga'+sep+'outputs'+sep+time+sep |
---|
| 18 | |
---|
| 19 | copy_code_files(output_dir,__file__) |
---|
| 20 | |
---|
[4926] | 21 | start_screen_catcher(output_dir,verbose=True) |
---|
| 22 | print'testing' |
---|
[4846] | 23 | |
---|
[4832] | 24 | data_dir='/d/ehn/4/ehn/tsunami/data/bathymetry/Cocos/bathymetry/' |
---|
| 25 | data_dir_name = data_dir + 'new_edited_cocos_final' |
---|
| 26 | |
---|
| 27 | data_dir_name_small=data_dir_name+'_small' |
---|
| 28 | data_dir_name_smaller=data_dir_name+'_smaller' |
---|
| 29 | data_dir_name_small_clip=data_dir_name+'_small_clip' |
---|
| 30 | data_dir_name_clip=data_dir_name+'_clip' |
---|
| 31 | #print'create Geospatial data2 objects from coast' |
---|
| 32 | #G = Geospatial_data(file_name = data_dir_name + '.txt') |
---|
| 33 | #print'start split' |
---|
| 34 | #G_small, G_other = G.split(factor=0.1, verbose=True) |
---|
| 35 | #print'start export' |
---|
| 36 | #G_small.export_points_file(data_dir_name_smaller + '.txt') |
---|
| 37 | |
---|
| 38 | mesh_file=data_dir+'cocos.msh' |
---|
| 39 | sww_file='test_final' |
---|
| 40 | |
---|
| 41 | #grid1_north_max=8783110.00 |
---|
| 42 | #grid1_north_min=8565314.80 |
---|
| 43 | #grid1_east_max=439850.00 |
---|
| 44 | #grid1_east_min=156494.00 |
---|
[4846] | 45 | grid1_res=370.4 |
---|
[4832] | 46 | # |
---|
| 47 | #grid2_north_max=8713178.04 |
---|
| 48 | #grid2_north_min=8625096.92 |
---|
| 49 | #grid2_east_max=294652.80 |
---|
| 50 | #grid2_east_min=224424.96 |
---|
[4846] | 51 | grid2_res=74.08 |
---|
[4832] | 52 | # |
---|
| 53 | #grid3_north_max=274295.91 |
---|
| 54 | #grid3_north_min=266767.44 |
---|
| 55 | #grid3_east_max=8665112.56 |
---|
| 56 | #grid3_east_min=8657346.91 |
---|
[4846] | 57 | grid3_res=14.82 |
---|
[4832] | 58 | # |
---|
| 59 | #grid4_north_max=270900.24 |
---|
| 60 | #grid4_north_min=269701.45 |
---|
| 61 | #grid4_east_max=8659120.05 |
---|
| 62 | #grid4_east_min=8659120.05 |
---|
[4846] | 63 | grid4_res=2.96 |
---|
[4832] | 64 | |
---|
| 65 | poly1 = [[439850.00,8565314.80],[439850.00, 8783110.00], |
---|
| 66 | [156494.00,8783110.00],[156494.00,8565314.80]] |
---|
| 67 | # 68450m2 |
---|
| 68 | |
---|
| 69 | poly2 = [[294652.80,8625096.92],[294652.80,8713178.04], |
---|
| 70 | [224424.96,8713178.04],[224424.96,8625096.92]] |
---|
| 71 | # 2738m2 |
---|
| 72 | |
---|
| 73 | poly3 = [[274295.91,8657346.91],[274295.91,8665112.56], |
---|
| 74 | [266767.44,8665112.56],[266767.44,8657346.91]] |
---|
| 75 | # 110m2 |
---|
| 76 | |
---|
| 77 | poly4 = [[270900.24,8659120.05],[270900.24,8660170.84], |
---|
| 78 | [269701.45,8660170.84],[269701.45,8659120.05]] |
---|
| 79 | # 4.3m2 |
---|
| 80 | |
---|
[4926] | 81 | res = 100 |
---|
[4832] | 82 | #G = Geospatial_data(file_name = data_dir_name + '.txt') |
---|
| 83 | #G1=G.clip(poly1) |
---|
| 84 | #G1.export_points_file(data_dir_name_clip + '.txt') |
---|
| 85 | |
---|
| 86 | |
---|
[4926] | 87 | #interior_polys = [[poly2,(grid2_res**2/2)*res],[poly3,(grid3_res**2/2)*res],[poly4,(grid4_res**2/2)*res]] |
---|
| 88 | interior_polys = [[poly2,(grid2_res**2)*2],[poly3,(grid3_res**2)*res],[poly4,(grid4_res**2)*res]] |
---|
[4832] | 89 | #interior_polys = [[poly4,4.3*res]] |
---|
| 90 | |
---|
| 91 | #create_mesh_from_regions(poly1, |
---|
| 92 | # boundary_tags={'side': [0,1,2,3]}, |
---|
| 93 | # maximum_triangle_area=68500*res, |
---|
| 94 | # interior_regions=interior_polys, |
---|
| 95 | # filename=mesh_file, |
---|
| 96 | # use_cache=True, |
---|
| 97 | # verbose=True) |
---|
[4926] | 98 | max_tri = (grid1_res**2/2)*res |
---|
| 99 | print 'max_tri',max_tri |
---|
[4832] | 100 | create_mesh_from_regions(poly1, |
---|
| 101 | #create_mesh_from_regions(poly3, |
---|
| 102 | boundary_tags={'side': [0,1,2,3]}, |
---|
[4926] | 103 | maximum_triangle_area=(grid1_res**2)*res, |
---|
| 104 | # maximum_triangle_area=(grid1_res**2)*1, |
---|
[4832] | 105 | interior_regions=interior_polys, |
---|
| 106 | filename=mesh_file, |
---|
[4993] | 107 | use_cache=False, |
---|
[4832] | 108 | verbose=True) |
---|
| 109 | |
---|
| 110 | |
---|
[4926] | 111 | #value, alpha = find_optimal_smoothing_parameter(data_file=data_dir_name_clip + '.txt', |
---|
| 112 | ## alpha_list=[0.0001,0.001,0.01,0.1,1.0], |
---|
| 113 | # alpha_list=[0.1], |
---|
| 114 | # boundary_poly=poly3, |
---|
| 115 | # mesh_file=mesh_file, |
---|
| 116 | # plot_name=None, |
---|
| 117 | # split_factor=0.1, |
---|
| 118 | # cache=True, |
---|
| 119 | # verbose=True) |
---|
[4832] | 120 | |
---|
[4926] | 121 | #print 'Setup computational domain', value, alpha |
---|
[4832] | 122 | |
---|
| 123 | #domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True) |
---|
| 124 | #above don't work |
---|
| 125 | domain = Domain(mesh_file, use_cache=True, verbose=True) |
---|
| 126 | |
---|
| 127 | print 'Start Set quantity' |
---|
| 128 | |
---|
| 129 | domain.set_quantity('elevation', |
---|
| 130 | filename = data_dir_name_clip + '.txt', |
---|
| 131 | use_cache = True, |
---|
| 132 | verbose = True, |
---|
[4926] | 133 | # alpha = alpha) |
---|
| 134 | alpha = 0.01) |
---|
[4832] | 135 | print 'Finished Set quantity' |
---|
| 136 | |
---|
| 137 | domain.set_name(sww_file) |
---|
[4846] | 138 | domain.set_datadir(output_dir) |
---|
[4926] | 139 | domain.set_default_order(1) # Apply second order scheme |
---|
[4832] | 140 | domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm |
---|
| 141 | domain.set_store_vertices_uniquely(False) |
---|
| 142 | domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) |
---|
| 143 | # domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK) |
---|
| 144 | domain.beta_h = 0 |
---|
| 145 | Br = Reflective_boundary(domain) |
---|
| 146 | |
---|
| 147 | |
---|
| 148 | domain.set_boundary({'side': Br}) |
---|
| 149 | |
---|
| 150 | for t in domain.evolve(yieldstep = 1,finaltime =1): |
---|
| 151 | domain.write_time() |
---|
| 152 | domain.write_boundary_statistics(tags = 'side') |
---|
| 153 | |
---|
[4846] | 154 | export_grid(output_dir+sww_file+".sww", extra_name_out = 'grid4', |
---|
[4832] | 155 | quantities = ['elevation'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation |
---|
| 156 | timestep = None, |
---|
| 157 | reduction = max, |
---|
[4926] | 158 | # cellsize = grid4_res*sqrt(1), |
---|
[4846] | 159 | cellsize = grid4_res*sqrt(res), |
---|
[4832] | 160 | NODATA_value = -9999, |
---|
| 161 | easting_min = 269701.45, |
---|
| 162 | easting_max = 270900.24, |
---|
| 163 | northing_min = 8659120.05, |
---|
| 164 | northing_max = 8660170.84, |
---|
| 165 | verbose = True, |
---|
| 166 | origin = None, |
---|
| 167 | datum = 'WGS84', |
---|
| 168 | format = 'asc') |
---|
| 169 | |
---|
[4846] | 170 | export_grid(output_dir+sww_file+".sww", extra_name_out = 'grid3', |
---|
[4832] | 171 | quantities = ['elevation'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation |
---|
| 172 | timestep = None, |
---|
| 173 | reduction = max, |
---|
[4926] | 174 | # cellsize = grid3_res*sqrt(1), |
---|
[4846] | 175 | cellsize = grid3_res*sqrt(res), |
---|
[4832] | 176 | NODATA_value = -9999, |
---|
| 177 | easting_min = 266767.44, |
---|
| 178 | easting_max = 274295.91, |
---|
| 179 | northing_min = 8657346.91, |
---|
| 180 | northing_max = 8665112.56, |
---|
| 181 | verbose = True, |
---|
| 182 | origin = None, |
---|
| 183 | datum = 'WGS84', |
---|
| 184 | format = 'asc') |
---|
| 185 | |
---|
[4846] | 186 | export_grid(output_dir+sww_file+".sww", extra_name_out = 'grid2', |
---|
[4832] | 187 | quantities = ['elevation'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation |
---|
| 188 | timestep = None, |
---|
| 189 | reduction = max, |
---|
[4926] | 190 | cellsize = grid2_res*sqrt(1), |
---|
| 191 | # cellsize = grid2_res*sqrt(res), |
---|
[4832] | 192 | NODATA_value = -9999, |
---|
| 193 | easting_min = 224424.96, |
---|
| 194 | easting_max = 294652.80, |
---|
| 195 | northing_min = 8625096.92, |
---|
| 196 | northing_max = 8713178.04, |
---|
| 197 | verbose = True, |
---|
| 198 | origin = None, |
---|
| 199 | datum = 'WGS84', |
---|
| 200 | format = 'asc') |
---|
| 201 | |
---|
[4846] | 202 | export_grid(output_dir+sww_file+".sww", extra_name_out = 'grid1', |
---|
[4832] | 203 | quantities = ['elevation'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation |
---|
| 204 | timestep = None, |
---|
| 205 | reduction = max, |
---|
[4926] | 206 | # cellsize = grid1_res*sqrt(1), |
---|
[4846] | 207 | cellsize = grid1_res*sqrt(res), |
---|
[4832] | 208 | NODATA_value = -9999, |
---|
| 209 | easting_min = 156494.00, |
---|
| 210 | easting_max = 439850.00, |
---|
| 211 | northing_min = 8565314.80, |
---|
| 212 | northing_max = 8783110.00, |
---|
| 213 | verbose = True, |
---|
| 214 | origin = None, |
---|
| 215 | datum = 'WGS84', |
---|
| 216 | format = 'asc') |
---|
| 217 | |
---|
| 218 | |
---|
| 219 | |
---|
| 220 | |
---|