"""Validation of the AnuGA implementation of the shallow water wave equation. This script sets up 2D circular island..... use 'res' parameter of 1 to make run fast. """ # Module imports from anuga.shallow_water.shallow_water_domain import Domain from anuga.shallow_water import Reflective_boundary from anuga.shallow_water import Transmissive_momentum_set_stage_boundary from anuga.shallow_water import Dirichlet_boundary, Time_boundary from anuga.utilities.file_utils import copy_code_files from anuga.shallow_water.sww_interrogate import get_maximum_inundation_data from anuga.abstract_2d_finite_volumes.util import file_function, \ sww2csv_gauges, csv2timeseries_graphs from anuga.pmesh.mesh_interface import create_mesh_from_regions from anuga.utilities.polygon import read_polygon from math import cos,pi,sin,tan import numpy as num from anuga.shallow_water.data_manager import load_csv_as_dict as csv2dict from time import localtime, strftime, gmtime from anuga.utilities.system_tools import get_user_name, get_host_name from os import sep, environ, getenv from anuga.config import netcdf_float from Scientific.IO.NetCDF import NetCDFFile #------------------------- # Create Domain from mesh #------------------------- user = get_user_name() home = getenv('INUNDATIONHOME') + sep +'data'+sep #Sandpit's parent dir time = strftime('%Y%m%d_%H%M%S',localtime()) res=0.05 dir_comment = time+'_res_'+str(res)+'_'+str(user) anuga_dir = home+'anuga_validation'+sep+'circular_island_tsunami_benchmark'+sep+'anuga'+sep output_dir = anuga_dir+'outputs'+sep+dir_comment+sep #output_dir = anuga_dir+'outputs'+sep+'20080222_062917_res_0.0025_nbartzis'+sep copy_code_files(output_dir,__file__) def get_xy(x=0,y=0,r=1,angle=1): #return x and y co-ordinates x1 = x + cos(angle)*r y1 = y + sin(-angle)*r # print "get",x,y,r,x1,y1,"angle",angle, cos(angle),sin(angle) return x1,y1, def prepare_timeboundary(filename): """Converting benchmark 2 time series to NetCDF tms file. This is a 'throw-away' code taylor made for files like 'Benchmark_2_input.txt' from the LWRU2 benchmark """ textversion = filename[:-4] + '.txt' print 'Preparing time boundary from %s' %textversion fid = open(textversion) #Skip first line line = fid.readline() #Read remaining lines lines = fid.readlines() fid.close() N = len(lines) T = num.zeros(N, num.float) #Time Q = num.zeros(N, num.float) #Values for i, line in enumerate(lines): fields = line.split() # print 'fields',i,line T[i] = float(fields[0]) Q[i] = float(fields[1]) #Create tms file print 'Writing to', filename fid = NetCDFFile(filename[:-4] + '.tms', 'w') fid.institution = 'Geoscience Australia' fid.description = 'Input wave for Benchmark 2' fid.starttime = 0.0 fid.createDimension('number_of_timesteps', len(T)) fid.createVariable('time', netcdf_float, ('number_of_timesteps',)) fid.variables['time'][:] = T fid.createVariable('stage', netcdf_float, ('number_of_timesteps',)) fid.variables['stage'][:] = Q[:] fid.createVariable('xmomentum', netcdf_float, ('number_of_timesteps',)) fid.variables['xmomentum'][:] = 0.0 fid.createVariable('ymomentum', netcdf_float, ('number_of_timesteps',)) fid.variables['ymomentum'][:] = 0.0 fid.close() # angle to make circular interior polygons angles1 = (0.0, 22.5, 45.0, 67.5, 90.0, 112.5, 135.0, 157.5, 180.0, 202.5, 225.0, 247.5, 270.0, 292.5, 315.0, 337.5) center_x = 12.96 center_y = 13.80 #L=3.6 and therefore the gauges 1,2,3 and 4 are located 1.8 from the toe of the island # 13.8 - 3.6(radii of island) - 1.8 = 8.4 L=8.4 #create polygons (circles) to have higher resolution poly = [] poly1 = [] internal_poly=[] for i,angle in enumerate(angles1): #convert Degs to Rads angle = ((angle-180)/180.0)*pi p = get_xy(center_x,center_y,4.5,angle) poly1.append([p[0],p[1]]) poly.append(poly1) internal_poly.append([poly1,.1*res]) poly2 = [] for i,angle in enumerate(angles1): #convert Degs to Rads angle = ((angle-180)/180.0)*pi p = get_xy(center_x,center_y,3.2,angle) poly2.append([p[0],p[1]]) poly.append(poly2) internal_poly.append([poly2,.01*res]) poly3 = [] for i,angle in enumerate(angles1): #convert Degs to Rads angle = ((angle-180)/180.0)*pi # p = get_xy(center_x,center_y,1.6,angle) p = get_xy(center_x,center_y,1.5,angle) poly3.append([p[0],p[1]]) poly.append(poly3) internal_poly.append([poly3,1*res]) #plot_polygons(poly,'poly.png') poly_all= [[0,L],[0,25],[30,25],[30,L]] create_mesh_from_regions(poly_all, boundary_tags={'wall': [0,2],'wave': [3], 'buffer':[1]}, maximum_triangle_area=1*res, interior_regions=internal_poly, filename='temp.msh', use_cache=False, verbose=True) domain = Domain('temp.msh', use_cache=True, verbose=True) print domain.statistics() #------------------------- # Initial Conditions #------------------------- #domain.set_quantity('friction', 0.0) domain.set_quantity('friction', 0.01)#for concrete #water_height = 0.322 water_height = 0.0 domain.set_quantity('stage', water_height) #attribute_dic, title_index_dic = csv2dict('sqrt_table.csv') def circular_island_elevation(x,y): water_depth = 0.32 list_xy = num.sqrt((center_x-x)**2+(center_y-L-y)**2) print 'x',min(x),max(x) print 'y',min(y),max(y) list_z=[] # for xy in list_xy_square: for distance in list_xy: # The cone has a 1/4 slope and a 7.2 diameter this make the height 0.9m # to get the stage at 0 elevation + cone height and - water depth z = -distance*.25+0.9-water_depth if z <0-water_depth: z=0-water_depth if z > .625-water_depth: z=0.625-water_depth list_z.append(z) return list_z domain.set_quantity('elevation', circular_island_elevation, verbose=True) ##------------------------- ## Set simulation parameters ##------------------------- model_name='wavetank_model' domain.set_name(model_name) # Name of output sww file domain.set_datadir(output_dir) domain.set_default_order(2) # Apply second order scheme #domain.set_all_limiters(0.9) # Max second order scheme (old lim) domain.set_minimum_storable_height(0.001) # Don't store w < 0.001m #domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK) domain.tight_slope_limiters = 1 sww_file=output_dir+sep+model_name+'.sww' #------------------------- # Boundary Conditions #------------------------- # Create boundary function from timeseries provided in file boundary_filename="ts2cnew1_input_20_80sec_new" prepare_timeboundary(boundary_filename+'.txt') function = file_function(boundary_filename+'.tms', domain, verbose=True) # Create and assign boundary objects Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function) Br = Reflective_boundary(domain) Bd = Dirichlet_boundary([water_height,0,0]) #domain.set_boundary({'wave': Bts, 'wall': Br, 'buffer':Bd}) domain.set_boundary({'wave': Bts, 'wall': Bd, 'buffer':Bd}) #domain.set_time(20.0) #------------------------- # Evolve through time #------------------------- import time t0 = time.time() for t in domain.evolve(yieldstep = 1, finaltime = 28): domain.write_time() # domain.write_time(track_speeds=False) for t in domain.evolve(yieldstep = 0.1, finaltime = 36 #for t in domain.evolve(yieldstep = 0.5, finaltime = 36 ,skip_initial_step = True): domain.write_time() for t in domain.evolve(yieldstep = 1, finaltime = 45 ,skip_initial_step = True): domain.write_time() print 'That took %.2f seconds' %(time.time()-t0) #define the segments to find the run-up height to compare to wave tank study angles = (0.0, 22.5, 45.0, 67.5, 75.0, 80.0, 85.0, 87.5, 90.0, 92.5, 95.0, 100.0, 105.0, 112.5, 135.0, 157.5, 180.0, 202.5, 225.0, 247.5, 270.0, 292.5, 315.0, 337.5) r1 = 1.5 r2 = 2.6 #degrees to check d=1 poly_segment = [] for i,angle in enumerate(angles): angle = ((angle-180)/180.0)*pi p1 = get_xy(center_x,center_y,r1,angle-0.01745*d) p2 = get_xy(center_x,center_y,r2,angle-0.01745*d) p3 = get_xy(center_x,center_y,r2,angle+0.01745*d) p4 = get_xy(center_x,center_y,r1,angle+0.01745*d) # print p1,p2,p3,p4,angle poly_segment =[[p1[0],p1[1]],[p2[0],p2[1]],[p3[0],p3[1]],[p4[0],p4[1]]] run_up, location = get_maximum_inundation_data(filename=sww_file,polygon=poly_segment, verbose=False) print 'maxd %.3f, %.6f, %.6f, %.5s, %.5s '%(((angle/pi)*180)+180, run_up*100, (run_up-water_height)*100, location[0],location[1]) #sww_file=anuga_dir+'outputs'+sep+'20080220_010628_res_0.05_nbartzis'+sep+'wavetank_model.sww' #sww_file=anuga_dir+'outputs'+sep+'20080219_233611_res_0.05_nbartzis'+sep+'wavetank_model.sww' sww2csv_gauges(sww_file=sww_file, gauge_file=anuga_dir+'gauges'+sep+'gauges.csv', # gauge_file='gauges.csv', quantities = ['stage','depth', 'elevation', 'xmomentum', 'ymomentum'], verbose=True, use_cache = True)