 Timestamp:
 Jun 15, 2010, 12:06:46 PM (13 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/anuga_core/source/anuga/shallow_water/test_data_manager.py
r7780 r7841 1121 1121 self.failUnless(0 ==1, 'Bad input did not throw exception error!') 1122 1122 1123 def test_export_grid_parallel(self):1124 """Test that sww information can be converted correctly to asc/prj1125 format readable by e.g. ArcView1126 """1127 1128 import time, os1129 from Scientific.IO.NetCDF import NetCDFFile1130 1131 base_name = 'tegp'1132 #Setup1133 self.domain.set_name(base_name+'_P0_8')1134 swwfile = self.domain.get_name() + '.sww'1135 1136 self.domain.set_datadir('.')1137 self.domain.format = 'sww'1138 self.domain.smooth = True1139 self.domain.set_quantity('elevation', lambda x,y: xy)1140 self.domain.set_quantity('stage', 1.0)1141 1142 self.domain.geo_reference = Geo_reference(56,308500,6189000)1143 1144 sww = SWW_file(self.domain)1145 sww.store_connectivity()1146 sww.store_timestep()1147 self.domain.evolve_to_end(finaltime = 0.0001)1148 #Setup1149 self.domain.set_name(base_name+'_P1_8')1150 swwfile2 = self.domain.get_name() + '.sww'1151 sww = SWW_file(self.domain)1152 sww.store_connectivity()1153 sww.store_timestep()1154 self.domain.evolve_to_end(finaltime = 0.0002)1155 sww.store_timestep()1156 1157 cellsize = 0.251158 #Check contents1159 #Get NetCDF1160 1161 fid = NetCDFFile(sww.filename, netcdf_mode_r)1162 1163 # Get the variables1164 x = fid.variables['x'][:]1165 y = fid.variables['y'][:]1166 z = fid.variables['elevation'][:]1167 time = fid.variables['time'][:]1168 stage = fid.variables['stage'][:]1169 1170 fid.close()1171 1172 #Export to ascii/prj files1173 extra_name_out = 'yeah'1174 sww2dem_batch(base_name,1175 quantities = ['elevation', 'depth'],1176 extra_name_out = extra_name_out,1177 cellsize = cellsize,1178 verbose = self.verbose,1179 format = 'asc')1180 1181 prjfile = base_name + '_P0_8_elevation_yeah.prj'1182 ascfile = base_name + '_P0_8_elevation_yeah.asc'1183 #Check asc file1184 ascid = open(ascfile)1185 lines = ascid.readlines()1186 ascid.close()1187 #Check grid values1188 for j in range(5):1189 L = lines[6+j].strip().split()1190 y = (4j) * cellsize1191 for i in range(5):1192 #print " i*cellsize  y", i*cellsize  y1193 #print "float(L[i])", float(L[i])1194 assert num.allclose(float(L[i]), i*cellsize  y)1195 #Cleanup1196 os.remove(prjfile)1197 os.remove(ascfile)1198 1199 prjfile = base_name + '_P1_8_elevation_yeah.prj'1200 ascfile = base_name + '_P1_8_elevation_yeah.asc'1201 #Check asc file1202 ascid = open(ascfile)1203 lines = ascid.readlines()1204 ascid.close()1205 #Check grid values1206 for j in range(5):1207 L = lines[6+j].strip().split()1208 y = (4j) * cellsize1209 for i in range(5):1210 #print " i*cellsize  y", i*cellsize  y1211 #print "float(L[i])", float(L[i])1212 assert num.allclose(float(L[i]), i*cellsize  y)1213 #Cleanup1214 os.remove(prjfile)1215 os.remove(ascfile)1216 os.remove(swwfile)1217 1218 #Check asc file1219 ascfile = base_name + '_P0_8_depth_yeah.asc'1220 prjfile = base_name + '_P0_8_depth_yeah.prj'1221 ascid = open(ascfile)1222 lines = ascid.readlines()1223 ascid.close()1224 #Check grid values1225 for j in range(5):1226 L = lines[6+j].strip().split()1227 y = (4j) * cellsize1228 for i in range(5):1229 assert num.allclose(float(L[i]), 1  (i*cellsize  y))1230 #Cleanup1231 os.remove(prjfile)1232 os.remove(ascfile)1233 1234 #Check asc file1235 ascfile = base_name + '_P1_8_depth_yeah.asc'1236 prjfile = base_name + '_P1_8_depth_yeah.prj'1237 ascid = open(ascfile)1238 lines = ascid.readlines()1239 ascid.close()1240 #Check grid values1241 for j in range(5):1242 L = lines[6+j].strip().split()1243 y = (4j) * cellsize1244 for i in range(5):1245 assert num.allclose(float(L[i]), 1  (i*cellsize  y))1246 #Cleanup1247 os.remove(prjfile)1248 os.remove(ascfile)1249 os.remove(swwfile2)1250 1251 1252 1253 def DISABLEDtest_sww2domain2(self):1254 ##################################################################1255 #Same as previous test, but this checks how NaNs are handled.1256 ##################################################################1257 1258 #FIXME: See ticket 2231259 1260 from mesh_factory import rectangular1261 1262 #Create basic mesh1263 points, vertices, boundary = rectangular(2,2)1264 1265 #Create shallow water domain1266 domain = Domain(points, vertices, boundary)1267 domain.smooth = False1268 domain.store = True1269 domain.set_name('test_file')1270 domain.set_datadir('.')1271 domain.default_order=21272 1273 domain.set_quantity('elevation', lambda x,y: x/3)1274 domain.set_quantity('friction', 0.1)1275 1276 from math import sin, pi1277 Br = Reflective_boundary(domain)1278 Bt = Transmissive_boundary(domain)1279 Bd = Dirichlet_boundary([0.2,0.,0.])1280 Bw = Time_boundary(domain=domain,1281 f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])1282 1283 domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})1284 1285 h = 0.051286 elevation = domain.quantities['elevation'].vertex_values1287 domain.set_quantity('stage', elevation + h)1288 1289 domain.check_integrity()1290 1291 for t in domain.evolve(yieldstep = 1, finaltime = 2.0):1292 pass1293 #domain.write_time()1294 1295 1296 filename = domain.datadir + os.sep + domain.get_name() + '.sww'1297 1298 # Fail because NaNs are present1299 #domain2 = sww2domain(filename,1300 # boundary,1301 # fail_if_NaN=True,1302 # verbose=self.verbose)1303 try:1304 domain2 = load_sww_as_domain(filename,1305 boundary,1306 fail_if_NaN=True,1307 verbose=self.verbose)1308 except DataDomainError:1309 # Now import it, filling NaNs to be 99991310 filler = 99991311 domain2 = load_sww_as_domain(filename,1312 None,1313 fail_if_NaN=False,1314 NaN_filler=filler,1315 verbose=self.verbose)1316 else:1317 raise Exception, 'should have failed'1318 1319 1320 # Now import it, filling NaNs to be 01321 filler = 99991322 domain2 = load_sww_as_domain(filename,1323 None,1324 fail_if_NaN=False,1325 NaN_filler=filler,1326 verbose=self.verbose)1327 1328 import sys; sys.exit()1329 1330 # Clean up1331 os.remove(filename)1332 1333 1334 bits = ['geo_reference.get_xllcorner()',1335 'geo_reference.get_yllcorner()',1336 'vertex_coordinates']1337 1338 for quantity in domain.quantities_to_be_stored:1339 bits.append('get_quantity("%s").get_integral()' %quantity)1340 bits.append('get_quantity("%s").get_values()' %quantity)1341 1342 for bit in bits:1343 # print 'testing that domain.'+bit+' has been restored'1344 assert num.allclose(eval('domain.'+bit),eval('domain2.'+bit))1345 1346 print1347 print1348 print domain2.get_quantity('xmomentum').get_values()1349 print1350 print domain2.get_quantity('stage').get_values()1351 print1352 1353 print 'filler', filler1354 print max(domain2.get_quantity('xmomentum').get_values().flat)1355 1356 assert max(max(domain2.get_quantity('xmomentum').get_values()))==filler1357 assert min(min(domain2.get_quantity('xmomentum').get_values()))==filler1358 assert max(max(domain2.get_quantity('ymomentum').get_values()))==filler1359 assert min(min(domain2.get_quantity('ymomentum').get_values()))==filler1360 1361 1362 1363 #FIXME This fails  smooth makes the comparism too hard for allclose1364 def ztest_sww2domain3(self):1365 ################################################1366 #DOMAIN.SMOOTH = TRUE !!!!!!!!!!!!!!!!!!!!!!!!!!1367 ################################################1368 from mesh_factory import rectangular1369 #Create basic mesh1370 1371 yiel=0.011372 points, vertices, boundary = rectangular(10,10)1373 1374 #Create shallow water domain1375 domain = Domain(points, vertices, boundary)1376 domain.geo_reference = Geo_reference(56,11,11)1377 domain.smooth = True1378 domain.store = True1379 domain.set_name('bedslope')1380 domain.default_order=21381 #Bedslope and friction1382 domain.set_quantity('elevation', lambda x,y: x/3)1383 domain.set_quantity('friction', 0.1)1384 # Boundary conditions1385 from math import sin, pi1386 Br = Reflective_boundary(domain)1387 Bt = Transmissive_boundary(domain)1388 Bd = Dirichlet_boundary([0.2,0.,0.])1389 Bw = Time_boundary(domain=domain,1390 f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])1391 1392 domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})1393 1394 domain.quantities_to_be_stored['xmomentum'] = 21395 domain.quantities_to_be_stored['ymomentum'] = 21396 #Initial condition1397 h = 0.051398 elevation = domain.quantities['elevation'].vertex_values1399 domain.set_quantity('stage', elevation + h)1400 1401 1402 domain.check_integrity()1403 #Evolution1404 for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):1405 # domain.write_time()1406 pass1407 1408 1409 filename = domain.datadir + os.sep + domain.get_name() + '.sww'1410 domain2 = load_sww_as_domain(filename,None,fail_if_NaN=False,verbose=self.verbose)1411 #points, vertices, boundary = rectangular(15,15)1412 #domain2.boundary = boundary1413 ###################1414 ##NOW TEST IT!!!1415 ###################1416 1417 os.remove(domain.get_name() + '.sww')1418 1419 #FIXME smooth domain so that they can be compared1420 1421 1422 bits = []1423 for quantity in domain.quantities_to_be_stored:1424 bits.append('quantities["%s"].get_integral()'%quantity)1425 1426 1427 for bit in bits:1428 #print 'testing that domain.'+bit+' has been restored'1429 #print bit1430 #print 'done'1431 #print ('domain.'+bit), eval('domain.'+bit)1432 #print ('domain2.'+bit), eval('domain2.'+bit)1433 assert num.allclose(eval('domain.'+bit),eval('domain2.'+bit),rtol=1.0e1,atol=1.e3)1434 pass1435 1436 ######################################1437 #Now evolve them both, just to be sure1438 ######################################x1439 domain.time = 0.1440 from time import sleep1441 1442 final = .51443 domain.set_quantity('friction', 0.1)1444 domain.store = False1445 domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Br})1446 1447 for t in domain.evolve(yieldstep = yiel, finaltime = final):1448 #domain.write_time()1449 pass1450 1451 domain2.smooth = True1452 domain2.store = False1453 domain2.default_order=21454 domain2.set_quantity('friction', 0.1)1455 #Bedslope and friction1456 # Boundary conditions1457 Bd2=Dirichlet_boundary([0.2,0.,0.])1458 Br2 = Reflective_boundary(domain2)1459 domain2.boundary = domain.boundary1460 #print 'domain2.boundary'1461 #print domain2.boundary1462 domain2.set_boundary({'left': Bd2, 'right': Bd2, 'top': Bd2, 'bottom': Br2})1463 #domain2.boundary = domain.boundary1464 #domain2.set_boundary({'exterior': Bd})1465 1466 domain2.check_integrity()1467 1468 for t in domain2.evolve(yieldstep = yiel, finaltime = final):1469 #domain2.write_time()1470 pass1471 1472 ###################1473 ##NOW TEST IT!!!1474 ##################1475 1476 print '><><><><>>'1477 bits = [ 'vertex_coordinates']1478 1479 for quantity in ['elevation','xmomentum','ymomentum']:1480 #bits.append('quantities["%s"].get_integral()'%quantity)1481 bits.append('get_quantity("%s").get_values()' %quantity)1482 1483 for bit in bits:1484 print bit1485 assert num.allclose(eval('domain.'+bit),eval('domain2.'+bit))1486 1123 1487 1124 def test_file_boundary_stsIV_sinewave_ordering(self):
Note: See TracChangeset
for help on using the changeset viewer.