source: anuga_validation/circular_island/run_circular.py @ 7746

Last change on this file since 7746 was 7746, checked in by hudson, 13 years ago

Upgraded validation to ANUGA 1.2

File size: 9.5 KB
RevLine 
[4777]1"""Validation of the AnuGA implementation of the shallow water wave equation.
2
3This script sets up 2D circular island.....
4
[5038]5use 'res' parameter of 1 to make run fast.
6
[4777]7"""
8
9# Module imports
10from anuga.shallow_water import Domain
11from anuga.shallow_water import Reflective_boundary
[7746]12from anuga.shallow_water import Transmissive_momentum_set_stage_boundary
[4777]13from anuga.shallow_water import Dirichlet_boundary, Time_boundary
[7746]14from anuga.utilities.file_utils import copy_code_files
15from anuga.shallow_water.data_manager import get_maximum_inundation_data
16from anuga.abstract_2d_finite_volumes.util import file_function, \
17                            sww2csv_gauges, csv2timeseries_graphs
[4777]18from anuga.pmesh.mesh_interface import create_mesh_from_regions
[7276]19from anuga.utilities.polygon import read_polygon
20from math import cos,pi,sin,tan
21import numpy as num
[7746]22from anuga.shallow_water.data_manager import load_csv_as_dict as csv2dict
[5065]23from time import localtime, strftime, gmtime
24from anuga.utilities.system_tools import get_user_name, get_host_name
25from os import sep, environ, getenv
[7276]26from anuga.config import netcdf_float
27from Scientific.IO.NetCDF import NetCDFFile
28
[4777]29#-------------------------
30# Create Domain from mesh
31#-------------------------
32
[5065]33user = get_user_name()
34home = getenv('INUNDATIONHOME') + sep +'data'+sep #Sandpit's parent dir 
35time = strftime('%Y%m%d_%H%M%S',localtime())
[5031]36
[5147]37res=0.05
[5065]38
39dir_comment = time+'_res_'+str(res)+'_'+str(user)
40
41anuga_dir = home+'anuga_validation'+sep+'circular_island_tsunami_benchmark'+sep+'anuga'+sep
42output_dir = anuga_dir+'outputs'+sep+dir_comment+sep
[5123]43#output_dir = anuga_dir+'outputs'+sep+'20080222_062917_res_0.0025_nbartzis'+sep
[5065]44copy_code_files(output_dir,__file__)
45
[4797]46def get_xy(x=0,y=0,r=1,angle=1):
[5031]47    #return x and y co-ordinates
[4797]48    x1 = x + cos(angle)*r
49    y1 = y + sin(-angle)*r
50#    print "get",x,y,r,x1,y1,"angle",angle, cos(angle),sin(angle)
51    return x1,y1,
[4777]52
53def prepare_timeboundary(filename):
54    """Converting benchmark 2 time series to NetCDF tms file.
55    This is a 'throw-away' code taylor made for files like
56    'Benchmark_2_input.txt' from the LWRU2 benchmark
57    """
58
59    textversion = filename[:-4] + '.txt'
60
61    print 'Preparing time boundary from %s' %textversion
62
63    fid = open(textversion)
64
65    #Skip first line
66    line = fid.readline()
67
68    #Read remaining lines
69    lines = fid.readlines()
70    fid.close()
71
72
73    N = len(lines)
[7276]74    T = num.zeros(N, num.float)  #Time
75    Q = num.zeros(N, num.float)  #Values
[4777]76
77    for i, line in enumerate(lines):
78        fields = line.split()
[4797]79#        print 'fields',i,line
[4777]80
81        T[i] = float(fields[0])
82        Q[i] = float(fields[1])
83
84
85    #Create tms file
86    print 'Writing to', filename
87    fid = NetCDFFile(filename[:-4] + '.tms', 'w')
88
89    fid.institution = 'Geoscience Australia'
90    fid.description = 'Input wave for Benchmark 2'
91    fid.starttime = 0.0
92    fid.createDimension('number_of_timesteps', len(T))
[7276]93    fid.createVariable('time', netcdf_float, ('number_of_timesteps',))
[4777]94    fid.variables['time'][:] = T
95
[7276]96    fid.createVariable('stage', netcdf_float, ('number_of_timesteps',))
[4777]97    fid.variables['stage'][:] = Q[:]
98
[7276]99    fid.createVariable('xmomentum', netcdf_float, ('number_of_timesteps',))
[4777]100    fid.variables['xmomentum'][:] = 0.0
101
[7276]102    fid.createVariable('ymomentum', netcdf_float, ('number_of_timesteps',))
[4777]103    fid.variables['ymomentum'][:] = 0.0
104
105    fid.close()
106
107
[5065]108# angle to make circular interior polygons
[4797]109angles1 = (0.0, 22.5, 45.0, 67.5, 90.0, 112.5, 135.0, 157.5, 180.0, 202.5, 225.0,
110          247.5, 270.0, 292.5, 315.0, 337.5)
[4777]111
[5031]112center_x = 12.96
113center_y = 13.80
114
[5038]115#L=3.6 and therefore the gauges 1,2,3 and 4 are located 1.8 from the toe of the island
116# 13.8 - 3.6(radii of island) - 1.8 = 8.4
117L=8.4
[4797]118
[4979]119#create polygons (circles) to have higher resolution
[4797]120poly = []
[4799]121poly1 = []
[4797]122internal_poly=[]
[4799]123for i,angle in enumerate(angles1):
124    #convert Degs to Rads
[4797]125    angle = ((angle-180)/180.0)*pi
126
[5147]127    p = get_xy(center_x,center_y,4.5,angle)
[4799]128    poly1.append([p[0],p[1]])
[4797]129
[4799]130poly.append(poly1)
[4877]131internal_poly.append([poly1,.1*res])
[4799]132
133poly2 = []
134for i,angle in enumerate(angles1):
135    #convert Degs to Rads
[4797]136    angle = ((angle-180)/180.0)*pi
[5147]137    p = get_xy(center_x,center_y,3.2,angle)
[4799]138    poly2.append([p[0],p[1]])
[4797]139   
[4799]140poly.append(poly2)
[4877]141internal_poly.append([poly2,.01*res])
[4797]142
[4799]143poly3 = []
144for i,angle in enumerate(angles1):
145    #convert Degs to Rads
[4797]146    angle = ((angle-180)/180.0)*pi
[5038]147#    p = get_xy(center_x,center_y,1.6,angle)
148    p = get_xy(center_x,center_y,1.5,angle)
[4799]149    poly3.append([p[0],p[1]])
[4797]150   
[4799]151poly.append(poly3)
[4877]152internal_poly.append([poly3,1*res])
[4797]153
[4987]154#plot_polygons(poly,'poly.png')
[4797]155
[5031]156poly_all= [[0,L],[0,25],[30,25],[30,L]]
[4777]157create_mesh_from_regions(poly_all,
[4797]158                         boundary_tags={'wall': [0,2],'wave': [3], 'buffer':[1]},
[5038]159                         maximum_triangle_area=1*res,
[4777]160                         interior_regions=internal_poly,
161                         filename='temp.msh',
162                         use_cache=False,
163                         verbose=True)
164
165domain = Domain('temp.msh', use_cache=True, verbose=True)
166print domain.statistics()
167
168#-------------------------
169# Initial Conditions
170#-------------------------
[5123]171#domain.set_quantity('friction', 0.0)
172domain.set_quantity('friction', 0.01)#for concrete
173#water_height = 0.322
174water_height = 0.0
175
[4777]176domain.set_quantity('stage', water_height)
177
[5031]178
[5123]179#attribute_dic, title_index_dic = csv2dict('sqrt_table.csv')
[4777]180
[5038]181def circular_island_elevation(x,y):
[5123]182    water_depth = 0.32
[7276]183    list_xy = num.sqrt((center_x-x)**2+(center_y-L-y)**2)
[5035]184    print 'x',min(x),max(x)
185    print 'y',min(y),max(y)
[5038]186    list_z=[]
[5123]187#    for xy in list_xy_square:
188    for distance in list_xy:
189
[5147]190#    The cone has a 1/4 slope and a 7.2 diameter this make the height 0.9m
191#    to get the stage at 0 elevation + cone height and - water depth
[5123]192        z = -distance*.25+0.9-water_depth
[5031]193
[5123]194        if z <0-water_depth:
195            z=0-water_depth
196        if z > .625-water_depth:
197            z=0.625-water_depth
[5038]198        list_z.append(z)
199    return list_z
[4980]200
[5038]201domain.set_quantity('elevation', circular_island_elevation, verbose=True)
[4980]202
203
[5123]204##-------------------------
205## Set simulation parameters
206##-------------------------
[5031]207model_name='wavetank_model'
208domain.set_name(model_name)  # Name of output sww file
[5065]209domain.set_datadir(output_dir)
[5123]210domain.set_default_order(2)               # Apply second order scheme
[4985]211#domain.set_all_limiters(0.9)              # Max second order scheme (old lim)
212domain.set_minimum_storable_height(0.001) # Don't store w < 0.001m
[5123]213#domain.set_maximum_allowed_speed(0.1)     # Allow a little runoff (0.1 is OK)
[4777]214domain.tight_slope_limiters = 1
[5123]215
[5065]216sww_file=output_dir+sep+model_name+'.sww'
[4777]217
218
219#-------------------------
220# Boundary Conditions
221#-------------------------
222
223# Create boundary function from timeseries provided in file
224
[5123]225boundary_filename="ts2cnew1_input_20_80sec_new"
[4777]226prepare_timeboundary(boundary_filename+'.txt')
227
228function = file_function(boundary_filename+'.tms',
229                         domain, verbose=True)
230
231# Create and assign boundary objects
232Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
233
234Br = Reflective_boundary(domain)
235Bd = Dirichlet_boundary([water_height,0,0])
[5123]236#domain.set_boundary({'wave': Bts, 'wall': Br, 'buffer':Bd})
237domain.set_boundary({'wave': Bts, 'wall': Bd, 'buffer':Bd})
[5038]238#domain.set_time(20.0)
[4777]239#-------------------------
240# Evolve through time
241#-------------------------
242import time
243t0 = time.time()
244
[5123]245for t in domain.evolve(yieldstep = 1, finaltime = 28):
[4777]246    domain.write_time()
247#    domain.write_time(track_speeds=False)
248
[4980]249for t in domain.evolve(yieldstep = 0.1, finaltime = 36
250#for t in domain.evolve(yieldstep = 0.5, finaltime = 36
[4797]251                       ,skip_initial_step = True): 
252    domain.write_time()
253   
[4985]254for t in domain.evolve(yieldstep = 1, finaltime = 45
[4797]255                       ,skip_initial_step = True): 
256    domain.write_time()
257
[4777]258print 'That took %.2f seconds' %(time.time()-t0)
[4797]259
260
261
[4979]262#define the segments to find the run-up height to compare to wave tank study
263angles = (0.0, 22.5, 45.0, 67.5, 75.0, 80.0, 85.0, 87.5, 90.0, 92.5,
264          95.0, 100.0, 105.0, 112.5, 135.0, 157.5, 180.0, 202.5, 225.0,
265          247.5, 270.0, 292.5, 315.0, 337.5)
[5038]266r1 = 1.5
267r2 = 2.6
268#degrees to check
[5065]269d=1
[4877]270
[4979]271poly_segment = []
272
[4797]273for i,angle in enumerate(angles):
[4979]274    angle = ((angle-180)/180.0)*pi
[5031]275    p1 = get_xy(center_x,center_y,r1,angle-0.01745*d)
276    p2 = get_xy(center_x,center_y,r2,angle-0.01745*d)
277    p3 = get_xy(center_x,center_y,r2,angle+0.01745*d)
278    p4 = get_xy(center_x,center_y,r1,angle+0.01745*d)
[4979]279#    print p1,p2,p3,p4,angle
280    poly_segment =[[p1[0],p1[1]],[p2[0],p2[1]],[p3[0],p3[1]],[p4[0],p4[1]]]
[4797]281   
[5065]282    run_up, location = get_maximum_inundation_data(filename=sww_file,polygon=poly_segment, verbose=False)
[4979]283   
[5038]284    print 'maxd %.3f, %.6f, %.6f, %.5s, %.5s '%(((angle/pi)*180)+180, run_up*100, (run_up-water_height)*100, location[0],location[1]) 
[4797]285
[5123]286#sww_file=anuga_dir+'outputs'+sep+'20080220_010628_res_0.05_nbartzis'+sep+'wavetank_model.sww'
287#sww_file=anuga_dir+'outputs'+sep+'20080219_233611_res_0.05_nbartzis'+sep+'wavetank_model.sww'
[5065]288sww2csv_gauges(sww_file=sww_file,
[5123]289                   gauge_file=anuga_dir+'gauges'+sep+'gauges.csv',
[5147]290#                   gauge_file='gauges.csv',
[5065]291                   quantities = ['stage','depth', 'elevation', 'xmomentum', 'ymomentum'],
292                   verbose=True,
293                   use_cache = True)
[4797]294
295
296
297
[4877]298
Note: See TracBrowser for help on using the repository browser.