"""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 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.shallow_water.data_manager import get_maximum_inundation_data, start_screen_catcher, copy_code_files 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#, plot_polygons from math import cos,pi,sin,tan,sqrt from Numeric import array, zeros, Float, allclose,resize from anuga.shallow_water.data_manager import 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 #------------------------- # 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 copy_code_files(output_dir,__file__) start_screen_catcher(output_dir) 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 from Numeric import array fid = open(textversion) #Skip first line line = fid.readline() #Read remaining lines lines = fid.readlines() fid.close() N = len(lines) T = zeros(N, Float) #Time Q = zeros(N, 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 from Scientific.IO.NetCDF import NetCDFFile 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', Float, ('number_of_timesteps',)) fid.variables['time'][:] = T fid.createVariable('stage', Float, ('number_of_timesteps',)) fid.variables['stage'][:] = Q[:] fid.createVariable('xmomentum', Float, ('number_of_timesteps',)) fid.variables['xmomentum'][:] = 0.0 fid.createVariable('ymomentum', 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 #d=0.75 #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.0,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,2.6,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) #water_height = 0.322 water_height = 0.0 domain.set_quantity('stage', water_height) attribute_dic, title_index_dic = csv2dict('sqrt_table.csv') #print attribute_dic #print attribute_dic['number'][10] def circular_island_elevation(x,y): list_xy_square = (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: # this finds the int that relates to the correct sqrt in the table nn=int(round(xy,2)*100) #the square root is to find the distance from the center to the current xy distance = float(attribute_dic['sqrt'][nn]) # zz=-float(attribute_dic['sqrt'][nn]) # print nn,zz # slope of cone is 1/4 and has a radii of 3.6 therefore the cone is 0.9 high. # So to have the elevation positive 0.90 has been added. Furthermore the #wavetank is submerged .32m so 0.9 - 0.32=.0578 # z = zz*.25+0.578 # z = -distance*.25+0.9 # z = -distance*.25+0.578 z = -distance*.25+0.9-0.322 # z = zz*.25+0.8975 if z <0-0.322: z=0-0.322 if z > .625-0.322: z=0.625-0.322 # print 'hello',z 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(1) # 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 domain.beta_h = 0.0 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]]] # print i, poly_segment, angle #print i,poly[i] 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='gauges.csv', quantities = ['stage','depth', 'elevation', 'xmomentum', 'ymomentum'], verbose=True, use_cache = True) #csv2timeseries_graphs(directories_dic={anuga_dir+sep+'boundaries'+sep:['wavetank',0, 0], # anuga_dir+'outputs'+sep+'20080219_233611_res_0.05_nbartzis'+sep:['Fixed Wave',0,0]}, # output_dir=anuga_dir+'outputs'+sep+'20080219_233611_res_0.05_nbartzis'+sep, # base_name='gauge_', # plot_numbers='', # quantities=['stage'], # extra_plot_name='', # assess_all_csv_files=True, # create_latex=False, # verbose=True)