Changeset 5031


Ignore:
Timestamp:
Feb 14, 2008, 1:50:57 PM (16 years ago)
Author:
nick
Message:

update circular island

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_validation/conical_island/run_circular.py

    r4987 r5031  
    1010from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
    1111from anuga.shallow_water import Dirichlet_boundary, Time_boundary
    12 from anuga.shallow_water.data_manager import get_maximum_inundation_data
     12from anuga.shallow_water.data_manager import get_maximum_inundation_data, start_screen_catcher
    1313from anuga.abstract_2d_finite_volumes.util import file_function, sww2csv_gauges
    1414from anuga.pmesh.mesh_interface import create_mesh_from_regions
     
    2323#-------------------------
    2424
     25#start_screen_catcher()
     26
    2527def get_xy(x=0,y=0,r=1,angle=1):
     28    #return x and y co-ordinates
    2629    x1 = x + cos(angle)*r
    2730    y1 = y + sin(-angle)*r
     
    9194          247.5, 270.0, 292.5, 315.0, 337.5)
    9295
    93 x = 12.96
    94 y = 13.80
     96#x = 12.96
     97#y = 13.80
     98center_x = 12.96
     99center_y = 13.80
     100
     101#wavelenght =
     102
    95103r1 = 1.5
    96104r2 = 2.5
     
    98106res=.05
    99107#res=1.
     108L=4.11
     109
    100110
    101111#create polygons (circles) to have higher resolution
     
    107117    angle = ((angle-180)/180.0)*pi
    108118
    109     p = get_xy(x,y,3.6,angle)
     119    p = get_xy(center_x,center_y,3.6,angle)
    110120    poly1.append([p[0],p[1]])
    111121
     
    117127    #convert Degs to Rads
    118128    angle = ((angle-180)/180.0)*pi
    119     p = get_xy(x,y,2.4,angle)
     129    p = get_xy(center_x,center_y,2.4,angle)
    120130    poly2.append([p[0],p[1]])
    121131   
     
    127137    #convert Degs to Rads
    128138    angle = ((angle-180)/180.0)*pi
    129     p = get_xy(x,y,1.6,angle)
     139    p = get_xy(center_x,center_y,1.6,angle)
    130140    poly3.append([p[0],p[1]])
    131141   
     
    137147    #convert Degs to Rads
    138148    angle = ((angle-180)/180.0)*pi
    139     p = get_xy(x,y,1.1,angle)
     149    p = get_xy(center_x,center_y,1.1,angle)
    140150    poly4.append([p[0],p[1]])
    141151   
     
    145155#plot_polygons(poly,'poly.png')
    146156
    147 poly_all= [[0,7],[0,25],[30,25],[30,7]]
     157#poly_all= [[0,7],[0,25],[30,25],[30,7]]
     158poly_all= [[0,L],[0,25],[30,25],[30,L]]
    148159create_mesh_from_regions(poly_all,
    149160                         boundary_tags={'wall': [0,2],'wave': [3], 'buffer':[1]},
     
    164175domain.set_quantity('stage', water_height)
    165176
    166 #attribute_dic, title_index_dic = csv2dict('sqrt_table.csv')
    167 attribute_dic, title_index_dic = csv2dict('sqrt_table_new.csv')
     177
     178attribute_dic, title_index_dic = csv2dict('sqrt_table.csv')
     179#print attribute_dic
    168180print attribute_dic['number'][10]
    169181
    170182def test(x,y):
    171     center_x = 12.96
     183#    center_x = 12.96
    172184#    center_y = 13.80
    173     center_y = 6.80
    174 
    175     du = (center_x-x)**2+(center_y-y)**2
     185#    center_y = 6.80
     186
     187    du = (center_x-x)**2+(center_y-L-y)**2
    176188#    print 'x',min(x),max(x)
    177189#    print 'y',min(y),max(y)
     
    182194        zz=-float(attribute_dic['sqrt'][nn])
    183195#        print nn,zz
    184         # this is the elevation for a cone with 1/4 slope and is ?? below the water .
     196
     197#    z = -square_root((center_x-x)**2+(center_y-y)**2)
    185198        z = zz*.25+0.8975
    186199
     
    189202        if z > .625:
    190203            z=0.625
    191         print 'hello',d,nn,zz,z
    192 #        print 'hello',d,n1,n2,n3,nn,zz,z
     204#        print 'hello',z
    193205        z1.append(z)
    194206    return z1
    195207
    196 domain.set_quantity('elevation',test, alpha=1., verbose=True)
     208domain.set_quantity('elevation',test, verbose=True)
    197209
    198210
     
    200212# Set simulation parameters
    201213#-------------------------
    202 domain.set_name('test1')  # Name of output sww file
     214model_name='wavetank_model'
     215domain.set_name(model_name)  # Name of output sww file
    203216domain.set_default_order(1)               # Apply second order scheme
    204217#domain.set_all_limiters(0.9)              # Max second order scheme (old lim)
     
    208221domain.beta_h = 0.0
    209222
    210 #Timings on AMD64-242 (beta_h=0)
    211 #  tight_slope_limiters = 0:
    212 #    3035s - 3110s
    213 #  tight_slope_limiters = 1:
    214 #    3000s - 3008s
    215 #
    216 # beta_h>0: In the order of 3200s
    217223
    218224#-------------------------
     
    223229#function = file_function(project.boundary_filename,
    224230#                         domain, verbose=True)
    225 def waveform(t):
    226         return (sin(t*2*pi))
    227231
    228232# Create and assign boundary objects
     
    236240Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
    237241
    238 #Bts = Transmissive_Momentum_Set_Stage_boundary(domain, waveform)
    239242Bw = Time_boundary(domain=domain,
    240243                   f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
     
    281284    angle = ((angle-180)/180.0)*pi
    282285#    angle = angle1*3.14157
    283     p1 = get_xy(x,y,r1,angle-0.01745*d)
    284     p2 = get_xy(x,y,r2,angle-0.01745*d)
    285     p3 = get_xy(x,y,r2,angle+0.01745*d)
    286     p4 = get_xy(x,y,r1,angle+0.01745*d)
     286    p1 = get_xy(center_x,center_y,r1,angle-0.01745*d)
     287    p2 = get_xy(center_x,center_y,r2,angle-0.01745*d)
     288    p3 = get_xy(center_x,center_y,r2,angle+0.01745*d)
     289    p4 = get_xy(center_x,center_y,r1,angle+0.01745*d)
    287290#    print p1,p2,p3,p4,angle
    288291    poly_segment =[[p1[0],p1[1]],[p2[0],p2[1]],[p3[0],p3[1]],[p4[0],p4[1]]]
     
    290293    #print i,poly[i]
    291294   
    292     run_up, x_y = get_maximum_inundation_data(filename='test.sww',polygon=poly_segment, verbose=False)
    293    
    294     print 'maximum_inundation_data',((angle/pi)*180)+180, run_up, x_y
    295 
    296 
    297 
    298 
    299 
    300 
    301 
     295    run_up, x_y = get_maximum_inundation_data(filename=model_name+'.sww',polygon=poly_segment, verbose=False)
     296   
     297    print 'maxd %.3f, %.6f, %.6f, %.5s, %.5s '%(((angle/pi)*180)+180, run_up*100, run_up-water_height*100, x_y[0],x_y[1])
     298
     299
     300
     301
     302
     303
     304
Note: See TracChangeset for help on using the changeset viewer.