1 | |
---|
2 | |
---|
3 | from anuga.geospatial_data.geospatial_data import * |
---|
4 | from os import sep, environ, getenv, getcwd,umask |
---|
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 |
---|
9 | from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files |
---|
10 | from math import sqrt |
---|
11 | from time import localtime, strftime, gmtime |
---|
12 | |
---|
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 | |
---|
21 | start_screen_catcher(output_dir,verbose=True) |
---|
22 | print'testing' |
---|
23 | |
---|
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 |
---|
45 | grid1_res=370.4 |
---|
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 |
---|
51 | grid2_res=74.08 |
---|
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 |
---|
57 | grid3_res=14.82 |
---|
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 |
---|
63 | grid4_res=2.96 |
---|
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 | |
---|
81 | res = 100 |
---|
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 | |
---|
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]] |
---|
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) |
---|
98 | max_tri = (grid1_res**2/2)*res |
---|
99 | print 'max_tri',max_tri |
---|
100 | create_mesh_from_regions(poly1, |
---|
101 | #create_mesh_from_regions(poly3, |
---|
102 | boundary_tags={'side': [0,1,2,3]}, |
---|
103 | maximum_triangle_area=(grid1_res**2)*res, |
---|
104 | # maximum_triangle_area=(grid1_res**2)*1, |
---|
105 | # maximum_triangle_area=110*res, |
---|
106 | interior_regions=interior_polys, |
---|
107 | filename=mesh_file, |
---|
108 | use_cache=True, |
---|
109 | verbose=True) |
---|
110 | |
---|
111 | |
---|
112 | #value, alpha = find_optimal_smoothing_parameter(data_file=data_dir_name_clip + '.txt', |
---|
113 | ## alpha_list=[0.0001,0.001,0.01,0.1,1.0], |
---|
114 | # alpha_list=[0.1], |
---|
115 | # boundary_poly=poly3, |
---|
116 | # mesh_file=mesh_file, |
---|
117 | # plot_name=None, |
---|
118 | # split_factor=0.1, |
---|
119 | # cache=True, |
---|
120 | # verbose=True) |
---|
121 | |
---|
122 | #print 'Setup computational domain', value, alpha |
---|
123 | |
---|
124 | #domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True) |
---|
125 | #above don't work |
---|
126 | domain = Domain(mesh_file, use_cache=True, verbose=True) |
---|
127 | |
---|
128 | print 'Start Set quantity' |
---|
129 | |
---|
130 | domain.set_quantity('elevation', |
---|
131 | filename = data_dir_name_clip + '.txt', |
---|
132 | use_cache = True, |
---|
133 | verbose = True, |
---|
134 | # alpha = alpha) |
---|
135 | alpha = 0.01) |
---|
136 | print 'Finished Set quantity' |
---|
137 | |
---|
138 | domain.set_name(sww_file) |
---|
139 | domain.set_datadir(output_dir) |
---|
140 | domain.set_default_order(1) # Apply second order scheme |
---|
141 | domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm |
---|
142 | domain.set_store_vertices_uniquely(False) |
---|
143 | domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) |
---|
144 | # domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK) |
---|
145 | domain.beta_h = 0 |
---|
146 | Br = Reflective_boundary(domain) |
---|
147 | |
---|
148 | |
---|
149 | domain.set_boundary({'side': Br}) |
---|
150 | |
---|
151 | for t in domain.evolve(yieldstep = 1,finaltime =1): |
---|
152 | domain.write_time() |
---|
153 | domain.write_boundary_statistics(tags = 'side') |
---|
154 | |
---|
155 | export_grid(output_dir+sww_file+".sww", extra_name_out = 'grid4', |
---|
156 | quantities = ['elevation'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation |
---|
157 | timestep = None, |
---|
158 | reduction = max, |
---|
159 | # cellsize = grid4_res*sqrt(1), |
---|
160 | cellsize = grid4_res*sqrt(res), |
---|
161 | NODATA_value = -9999, |
---|
162 | easting_min = 269701.45, |
---|
163 | easting_max = 270900.24, |
---|
164 | northing_min = 8659120.05, |
---|
165 | northing_max = 8660170.84, |
---|
166 | verbose = True, |
---|
167 | origin = None, |
---|
168 | datum = 'WGS84', |
---|
169 | format = 'asc') |
---|
170 | |
---|
171 | export_grid(output_dir+sww_file+".sww", extra_name_out = 'grid3', |
---|
172 | quantities = ['elevation'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation |
---|
173 | timestep = None, |
---|
174 | reduction = max, |
---|
175 | # cellsize = grid3_res*sqrt(1), |
---|
176 | cellsize = grid3_res*sqrt(res), |
---|
177 | NODATA_value = -9999, |
---|
178 | easting_min = 266767.44, |
---|
179 | easting_max = 274295.91, |
---|
180 | northing_min = 8657346.91, |
---|
181 | northing_max = 8665112.56, |
---|
182 | verbose = True, |
---|
183 | origin = None, |
---|
184 | datum = 'WGS84', |
---|
185 | format = 'asc') |
---|
186 | |
---|
187 | export_grid(output_dir+sww_file+".sww", extra_name_out = 'grid2', |
---|
188 | quantities = ['elevation'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation |
---|
189 | timestep = None, |
---|
190 | reduction = max, |
---|
191 | cellsize = grid2_res*sqrt(1), |
---|
192 | # cellsize = grid2_res*sqrt(res), |
---|
193 | NODATA_value = -9999, |
---|
194 | easting_min = 224424.96, |
---|
195 | easting_max = 294652.80, |
---|
196 | northing_min = 8625096.92, |
---|
197 | northing_max = 8713178.04, |
---|
198 | verbose = True, |
---|
199 | origin = None, |
---|
200 | datum = 'WGS84', |
---|
201 | format = 'asc') |
---|
202 | |
---|
203 | export_grid(output_dir+sww_file+".sww", extra_name_out = 'grid1', |
---|
204 | quantities = ['elevation'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation |
---|
205 | timestep = None, |
---|
206 | reduction = max, |
---|
207 | # cellsize = grid1_res*sqrt(1), |
---|
208 | cellsize = grid1_res*sqrt(res), |
---|
209 | NODATA_value = -9999, |
---|
210 | easting_min = 156494.00, |
---|
211 | easting_max = 439850.00, |
---|
212 | northing_min = 8565314.80, |
---|
213 | northing_max = 8783110.00, |
---|
214 | verbose = True, |
---|
215 | origin = None, |
---|
216 | datum = 'WGS84', |
---|
217 | format = 'asc') |
---|
218 | |
---|
219 | |
---|
220 | |
---|
221 | |
---|